i Seismotectonic and seismic hazard studies for Madagascar by Tsiriandrimanana Rakotondraibe A thesis submitted to the Faculty of Science, University of the Witwatersrand, Johannesburg, in fulfillment of the requirements for the degree of Doctor of Philosophy ii iii Declaration I declare that this thesis is my individual work. It is submitted in fulfillment of the requirements for the Degree of Doctor of Philosophy in the University of the Witwatersrand, Johannesburg. This thesis was not submitted before for any degree or examination in any other University. Date: 24th October 2022 iv v Abstract Madagascar is an island that was once situated in the centre of the Gondwana supercontinent but is now located 1000 km from the East African Rift and about 2000 km from the mid-ocean ridge of the Indian Ocean. Its separation from the African continent (ca. 165-130 Ma) and the Indian sub-continent (ca. 90 Ma) and mantle upwelling that gave rise to the recent volcanism (ca. 28-0.5 Ma) has had significant impact on the geological structure and neotectonics of the island. This thesis focuses on the seismotectonic evolution of Madagascar and assesses the seismic hazard. Data for this project was derived from the 28-station MAdagascar – COmoros-MOzambique (MACOMO) array, as well as stations from the concurrent SEismological signatures in the Lithosphere/Asthenosphere system of SOuthern MAdagascar (SELASOMA) and Réunion Hotspot and Upper Mantle - Réunions Unterer Mantel (RHUM-RUM) projects, and permanent stations of the national seismological network. The combined 56 stations provided a better distribution and density than any previous study. The thesis is presented in three main parts. The first part focuses on seismicity and focal mechanisms. The 695 events recorded and located during the 23-month period of MACOMO observations in 2011-2013 have magnitudes that vary from M1 to M5.3. Twenty-three focal mechanism solutions were determined using P- wave first motion polarities and some clear S-wave polarities. In agreement with previous studies, it was found that central Madagascar is the most seismically- active region. Most earthquakes are located near known faults. Normal faulting dominates the central area, while strike-slip and reverse fault mechanisms are predominantly observed near the Antsaba and Sandrakota shear zones in the north, and the Bekily and Ampanihy shear zones in the south. The focal mechanism studies revealed extensional faulting in central Madagascar, which can be explained by the reactivation of older faults by a thermal bulge underlying the region. The faulting style was not well defined in the northern and southern parts as few solutions were determined due to the lower density of stations. vi In the second part of this thesis, 182 focal mechanism solutions were used to study the characteristics of the stress field. Calculations were performed using a linear inversion technique based on the Michael (1984) method and using a recently calculated velocity model of Madagascar. Results give shape ratios (i.e., relative magnitude of stress axes) that vary from 0.7 to 0.9; a maximum principal stress axis (σ1) that is near-vertical and a minimum principal stress axis (σ3) that is near- horizontal with a NW-SE orientation. This corresponds to a normal faulting regime. The regional variation of stress and the characteristics of the shape ratio and friction values were determined. The average direction of the stress field is approximately NW-SE in central Madagascar, and trends W-E in the southern and northern parts of the island. In central Madagascar the shape ratio and friction vary from 0.83 to 0.97 and 0.45 to 0.65, respectively; while in the northern and southern parts of Madagascar they vary from 0.3 to 0.5 and 0.5 to 0.7, respectively. The third part of the thesis reports on the assessment of the seismic hazard in Madagascar using data from a catalogue that combines data from the NDC (National Data Center), bulletins produced by the current study for the period between 1975 and 2016, as well as historical earthquake data. The analysis reveals that Madagascar has moderate seismic activity with local magnitude ML ≤ 6. The relationship between the moment magnitude MW and local magnitude ML values was estimated in order to unify magnitude characteristics and create the homogeneous catalogue that is required for seismic hazard studies. A Probabilistic Seismic Hazard Assessment approach (PSHA) was performed using the Kijko and Sellevoll (KS) and Cornell-McGuire methods with two different ground motion models. In order to implement the alternative parameter values in seismic hazard calculation and to attribute uncertainties associated to the input parameter, a logic tree formalism was applied. The Peak Ground Acceleration (PGA) values near major towns in the central region for 10% and 2% probability of occurrence in 50 years are 0.04 g and 0.11 g, respectively. The spectral acceleration for 1.0 s and 3.0 s vary from 0.001 g to 0.01 g and 0.001 g to 0.005 g for 10% probability of vii occurrence in 50 years, respectively; and from 0.02 g to 0.10 g and 0.002 g to 0.01 g for 2% probability of occurrence in 50 years. This work demonstrates that the deployment of additional temporary seismic stations significantly improved our understanding of the seismotectonics of Madagascar, particularly in the northern and southern regions. In summary, we obtained a new magnitude relation, more accurate epicenter locations, better focal mechanism descriptions, and improved the assessment of seismic hazard. viii Dedication In memory of my mother Rasoanirina Justine Claire 1952-2012 ix x Acknowledgements Primary, I would like to thank God for giving me strength during this research. This PhD thesis has been done mostly at the School of Geosciences, University of the Witwatersrand Johannesburg, South Africa under guidance of Prof. Raymond Durrheim and Prof. Andrew Nyblade. It was a great work under their consideration; I really appreciate their concern, comment, and advice. I would like to give special thanks to Prof. Michael Wysession for his help, encouragement, and acceptance of me as member on his team during the deployment of the MACOMO stations. I would like to thank Prof. Gerard Rambolamanana for giving me the opportunity to participate to the deployment of the MACOMO stations and for his encouragements and advice. I would like to thank also Profs Frederic Tillman, Karin Sigloch and Guilhem Barruol for allowing me to use the data from the SELASOMA and RHUM-RUM projects. I am so grateful to all of the field team members during the MAdagascar- COmoros-Mozambique (MACOMO) installation, Ghassan Aleqabi, Patrick Shore, Martin Pratt, Fenitra Andriampenomanana, Rina Randrianasolo and the drivers. I am so grateful to Mrs. Sharon Ellis, Mrs. Antonia Money and Mrs. Barbara Peragine, the administrative assistants, for their support and assistance during my study period. I acknowledge all the members of the AfricaArray group at the Department of Geosciences, School of Earth and Mineral Sciences, Penn State University, USA. I am grateful to all members of the group meeting in the Seismic Lab at the School of Geosciences, University of the Witwatersrand, Johannesburg, South Africa, for their support and encouragement. I would like to thank the School of Geosciences as a whole; the staff including the head of the Department of Geosciences. xi I would like to thank Michele Chernega for her support during my stay in State College while at PennState University, USA. I am heartily thankful for her encouragement, advice and for care and everything and my friends Carol Destro, Atefeh Mohammadpour, Laiz Fiorilli, Jocely Souza for making it memorable. I would like to thank all the staff at the Institut et Observatoire de Géophysique d'Antananarivo (IOGA), all the members at the "Laboratoire de Sismologie et Infrason (LSI)" at the IOGA and the University of Antananarivo, Madagascar for their encouragement and help during the preparation of this study. I would like to thank Dr Andriamiranto Raveloson, Dr Elisa Rindraharisaona, Dr Fenitra Andriampenomanana and Khile Homann for their help with scripts and advice. I would like to thank all of my friends from Madagascar, South Africa and the USA that helped me to accomplish this study. I would like to thank all of my family members, especially for their continuous support and encouragement over the years. Their support, love and encouragement have helped me during the challenging times, and also my mom who passed away before I started this study. This study was funded by NSF awards EAR-0838426 and 0838387, and supported by the AfricaArray Program. Prof. Raymond Durrheim acknowledges the support of the National Research Foundation through the South African Research Chairs Initiative. xii List of publications Peer-reviewed publication Andriampenomanana, F., Nyblade, A.A., Wysession, M.E., Durrheim, R.J., Tilmann, F., Barruol, G., Rambolamanana, G., and Rakotondraibe, T. (2020). Seismic velocity and anisotropy of the uppermost mantle beneath Madagascar from Pn tomography, Geophysical Journal International, https://doi.org/10.1093/gji/ggaa458 Rakotondraibe T., Nyblade,A.A., Wysession, M.E., Durrheim, R.J., Rambolamanana, G., Aleqabi, G., Shore, P.J., Pratt, M.J., Andriampenomanana, F., Rümpker, G., and Rindraharisaona, E. (2020). Seismicity and seismotectonics of Madagascar revealed by the 2011- 2013 deployment of the island-wide MACOMO broadband seismic array, Tectonophysics, 790, 228547, https://doi.org/10.1016/j.tecto.2020.228547 Ramirez, C., Nyblade, A.A., Wysession, M.E., Pratt, M.J., Andriampenomanana, F., and Rakotondraibe, T. (2018). Complex seismic anisotropy in Madagascar revealed by shear wave splitting measurements, Geophysical Journal International, 215(3), 1718-1727, https://doi.org/10.1093/gji/ggy367 Andriampenomanana, F., Nyblade, A.A., Wysession, M.E., Durrheim, R.J., Tilmann, F., Julià, J., Pratt, M.J., Rambolamanana, G., Aleqabi, G., Shore, P.J., and Rakotondraibe T. (2017). The structure of the crust and uppermost mantle beneath Madagascar, Geophysical Journal International, 210 (3), 1525 - 1544,http://doi:10.1093/gji/ggx243 Pratt, M.J., Wysession, M.E., Aleqabi, G.I., Wiens, D.A., Nyblade, A.A, Shore, P. J., Rambolamanana, G., Andriampenomanana, F., Rakotondraibe, T., Tucker, R.D., Barruol, G., and Rindraharisaona, E. (2017). Shear-velocity structure of the crust and uppermantle of Madagascar derived from surface wave tomography, Earth and Planetary Science Letters, 458, 405 - 417, http://doi:10.1016/j.epsl.2016.10.041 Peer-reviewed conference proceedings 06 - 09 October 2019: 16th Biennial Conference and Exhibition of the South African Geophysical Association (SAGA) (South Africa). Oral presentation: “Seismic hazard studies in Madagascar”. https://doi.org/10.1093/gji/ggaa458 https://doi.org/10.1016/j.tecto.2020.228547 https://doi.org/10.1093/gji/ggy367 http://doi:10.1016/j.epsl.2016.10.041 xiii 07 - 12 July 2019: 14th Annual AfricaArray Scientific Workshop at the School of Geosciences of the University of Witwatersrand (South Africa). Oral presentation: “Seismic hazard studies for Madagascar”. 02 – 07 September 2018: 36th General Assembly of the European Seismological Commission in Valletta, Malta (ESC2018-S20-435). Oral and Poster presentation: “Preliminary results of stress derivation from earthquake focal mechanisms in Madagascar”. 25 June – 26 June 2018: 13th Annual Africa Array Scientific Workshop at the School of Geosciences of the University of Witwatersrand (South Africa). Poster Presentation “Preliminary results of stress derivation from earthquake focal mechanisms in Madagascar”. 11 - 15 December 2017: Fall Meeting of the American Geophysical Union, New Orleans, USA. Wysession, ME, Andriampenomanana, F, Nyblade, AA, Durrheim, R, Wiens, DA, Rambolamanana, G, Pratt, MJ, Aleqabi, G, Shore, PJ, Rakotondraibe, T, Tucker, RD and Tilmann, F, 2017, T41E-03: Normal- Faulting in Madagascar: Another Round of Continental Rifting? 13 - 14 November 2017: National Data Center (NDC) Capacity building Workshop and RSTT model (Namibia). Poster presentation: “Seismotectonic studies in Madagascar”. 10 -13 September 2017: 15th Biennial Conference and Exhibition of the South African Geophysical Association (SAGA) (South Africa). Oral presentation: “Focal mechanism and stress field inversion in the central part of Madagascar”. 29 June – 30 June 2017: 12th Annual AfricaArray Scientific Workshop at the School of Geosciences of the University of Witwatersrand (South Africa). Oral presentation: “Probabilistic seismic hazard for Madagascar”. 27 August - 04 September 2016: 35th International Geological Congress (IGC) (South Africa). Poster presentation: “Seimotectonic studies for Madagascar”. 02-05 April 2016: 1 st General Assembly of the African Seismological Commission (AfSC), (Egypt). Poster presentation: “Seismotectonics studies for Madagascar”. 17 - 23 January 2016: 11th Annual AfricaArray Scientific Workshop at the School of Geosciences of the University of Witwatersrand (South Africa). Oral presentation: “Determination of Fault structure for Madagascar”. 06 - 09 September 2015: 14th Biennial Conference and Exhibition of the South African Geophysical Association (SAGA) (South Africa). Oral presentation: “Seismicity and Local Magnitude scale for Madagascar”. Champagne Castle, SOUTH AFRICA, pp. 173-176. xiv 19 - 20 January 2015: 10th Annual AfricaArray Scientific Workshop at the School of Geosciences of the University of Witwatersrand (South Africa). Oral presentation: “Seismicity map of Madagascar using data from MACOMO project”. xv Contents Declaration ................................................................................................................... iii Abstract ......................................................................................................................... v Dedication................................................................................................................... viii Acknowledgements ....................................................................................................... x List of publications ..................................................................................................... xii List of figures .............................................................................................................. xx List of tables ............................................................................................................. xxvi List of Acronyms ................................................................................................... xxviii Chapter 1: General Introduction ................................................................................ 1 1.1 Introduction.................................................................................................... 1 1.2 Objectives ...................................................................................................... 2 1.3 Outline ........................................................................................................... 2 1.3.1 Seismicity studies .................................................................................... 2 1.3.2 Fault mechanisms and stress field model ................................................ 3 1.3.3 Magnitude scale and seismic hazard assessment .................................... 3 1.4 Problem statement ......................................................................................... 3 1.5 Literature review: geological setting ............................................................. 4 1.5.1 Precambrian basement ............................................................................ 6 1.5.2 Sedimentary cover................................................................................... 9 1.5.3 Volcanic fields ........................................................................................ 9 1.6 Literature review: previous geophysical work ............................................ 11 1.6.1 Seismicity studies of Madagascar ......................................................... 12 1.6.2 Magnitude scale .................................................................................... 14 1.6.3 Velocity structure .................................................................................. 15 1.6.4 Focal mechanisms ................................................................................. 18 1.6.5 Hazard studies ....................................................................................... 19 Chapter 2: Theory ...................................................................................................... 20 xvi 2.1 Seismic waves .............................................................................................. 20 2.2 Source parameters and mechanisms ............................................................ 22 2.2.1 Location and origin time ....................................................................... 22 2.2.2 Correction for instrument response ....................................................... 26 2.2.3 Fault geometry ...................................................................................... 27 2.2.4 Focal mechanism solutions ................................................................... 29 2.3 Stress field ................................................................................................... 37 2.3.1 Basic concepts ....................................................................................... 37 2.3.2 Mohr’s circle of stress ........................................................................... 38 2.3.3 P- and T-axis ......................................................................................... 41 2.3.4 Stress inversion ..................................................................................... 41 2.4 Earthquake magnitudes ................................................................................ 47 2.4.1 Local magnitude .................................................................................... 47 2.4.2 Duration magnitude ............................................................................... 50 2.4.3 Moment magnitude ............................................................................... 51 2.4.4 Relation between MW, ML, MS, mB and mb ............................................ 54 2.5 Seismic Hazard Assessment ........................................................................ 55 2.5.1 Basic concepts ....................................................................................... 55 2.5.2 Cornell-McGuire approach ................................................................... 58 2.5.3 Kijko and Sellevoll approach (KS) ....................................................... 62 2.5.4 Logic tree .............................................................................................. 68 2.6 Computational packages .............................................................................. 69 2.6.1 ANTELOPE software ........................................................................... 69 2.6.2 HYPOELLIPSE .................................................................................... 70 2.6.3 HYPODD .............................................................................................. 72 2.6.5 STRESSINVERSE program ................................................................. 74 2.6.6 HA3 program ........................................................................................ 75 2.6.7 CRISIS software ................................................................................... 75 Chapter 3: Data and Methodology ........................................................................... 77 3.1 Data .............................................................................................................. 77 3.1.2 Permanent seismic network .................................................................. 77 xvii 3.1.3 Temporary seismic networks ................................................................ 79 3.2 Methodology ................................................................................................ 90 3.2.1 Completeness of the catalogue .............................................................. 90 3.2.2 Declustering method ............................................................................. 90 3.2.3 Seismic source zone .............................................................................. 92 3.2.4 Earthquake catalogue ............................................................................ 92 3.2.5 Magnitude homogenization ................................................................... 93 Chapter 4: Seismotectonics of Madagascar ............................................................. 98 4.1 Introduction.................................................................................................. 98 4.2 Seismicity .................................................................................................... 99 4.2.1 Location of earthquakes ........................................................................ 99 4.2.2 Northern region ................................................................................... 102 4.2.3 Central region ...................................................................................... 103 4.2.4 Southern region ................................................................................... 107 4.3 Focal mechanisms ...................................................................................... 109 4.3.1 Computation ........................................................................................ 109 4.3.2 Discussion ........................................................................................... 113 4.4 Magnitude calibration ................................................................................ 117 4.5 Discussion .................................................................................................. 120 4.6 Conclusion ................................................................................................. 123 Chapter 5: Stress field inversion ............................................................................. 125 5.1 Introduction................................................................................................ 125 5.2 Focal mechanisms ...................................................................................... 125 5. 3 Stress Inversion ......................................................................................... 130 5.4 Results of inversion ................................................................................... 131 5.4.1 Bemarivo ............................................................................................. 132 5.4.2 Bekily .................................................................................................. 133 5.4.3 Famoizankova ..................................................................................... 135 5.4.4 Alaotra ................................................................................................. 136 5.4.5 Ankaratra ............................................................................................. 141 5.4.6 Itremo .................................................................................................. 145 xviii 5.5 Discussion .................................................................................................. 148 5.5.1 Northern region ................................................................................... 148 5.5.2 Central region ...................................................................................... 149 5.5.3 Bekily region ....................................................................................... 150 5.6 Conclusion ................................................................................................. 151 Chapter 6: Seismic hazard analysis ........................................................................ 152 6.1 Introduction................................................................................................ 152 6.2 Methodology .............................................................................................. 153 6.2.1 Declustering ........................................................................................ 158 6.2.2 Definition and characterization of source zones ................................. 165 6.2.3 Ground Motion Prediction Equations ................................................. 174 6.3 Results ....................................................................................................... 181 6.3.1 Hazard maps in term of PGA .............................................................. 181 6.3.2 Uniform hazard spectrum .................................................................... 184 6.4 Discussions and conclusion ....................................................................... 185 6.4.1 Northern and southern part of Madagascar ......................................... 187 6.4.2 Central part of Madagascar ................................................................. 188 Chapter 7: Discussion and Conclusion ................................................................... 189 7.1 Introduction................................................................................................ 189 7.2 General discussion and conclusion ............................................................ 190 7.2.1 Seismicity and magnitude ................................................................... 190 7.2 2 Focal mechanisms ............................................................................... 191 7.2.3 Stress field ........................................................................................... 193 7.2.4 Seismic hazard .................................................................................... 195 7. 3 Recommendations for future studies ........................................................ 197 References ................................................................................................................. 199 Appendix A: Seismotectonic studies of Madagascar (Published Paper) ............. 220 Appendix B: Example of input files for HYPOELLIPSE and HYPODD .......... 245 Appendix C: Earthquake location and relocation lists ......................................... 247 Appendix D: Example of FOCMEC input files ..................................................... 273 xix Appendix E: Paper about Stress field in Madagascar paper (In Prep.) ............. 277 Appendix F. Focal mechanism list for stress determination ................................ 294 Appendix G: Paper about Seismic hazard (In prep.) ........................................... 300 Appendix H: Moment magnitude determination .................................................. 315 xx List of figures Figure 1. 1: Gondwana showing the position of Madagascar at the end of the Neoproterozoic Era (~550 Ma) (Kusky et al., 2003). ..................................................... 5 Figure 1. 2: Displacement of Madagascar along the Davie Fracture Zone (Mascle et al., 1987). ..................................................................................................................... 5 Figure 1. 4: Ankaratra and Lac Itasy volcanic fields in central Madagascar (after Bésairie, 1964). ............................................................................................................. 11 Figure 1. 5: Seismicity and focal mechanisms determined by earlier studies.Red circles are events in the ISC Bulletin with magnitudes M>4. ....................................... 13 Figure 1. 6: (a) S-wave velocity model from all stations (Andriampenomanana et al. 2017), (b) Calculated vertical P-wave velocity structure for Madagascar using results from the (c) Wadati plot. ................................................................................... 16 Figure 2. 1: Main types of waves produced by an earthquake. (https://courses.lumenlearning.com/geo/chapter/reading-studying-the-earths- interior ) ......................................................................................................................... 21 Figure 2. 2: Vertical component seismogram recorded at the SBV seismic station produced by a M3.5 earthquake on March 16th, 2012 at 11:56:55 UTC. ..................... 21 Figure 2. 3: Earthquake focus (hypocenter) and epicenter. .......................................... 22 Figure 2. 4: Triangulation method to locate earthquakes. http://www.bbc.co.uk/education/guides/zq7ycdm/revision/2 )..................................... 23 Figure 2. 5: Geometrical description of a fault plane 𝜸 is the rake direction and n is the normal to the fault plane.......................................................................................... 28 Figure 2. 6: Illustration of the method for calculating strike, dip, and rake angles ...... 29 Figure 2. 7: Measurement of the azimuth angle ............................................................ 30 Figure 2. 8: Definition of the take-off angle ................................................................. 30 Figure 2. 9: Illustration of P-wave upward and downward first motions (Havskov and Ottemöller, 2010) ................................................................................................... 31 Figure 2. 10: Transformation of ZNE to ZRT and LQT ............................................... 32 xxi Figure 2. 11: Examples of polarity picking on a local M3.1 event recorded at station BKTA on 4 July 2013 at 13:40:59.08. .............................................................. 33 Figure 2. 12: Fault mechanisms and ‘beach ball’ diagrams. ......................................... 34 Figure 2. 13: Focal mechanism for a normal faulting earthquake with the actual fault plane and auxiliary plane in a section BB'. ........................................................... 35 Figure 2. 14: Normal fault plane solution showing the P- and T-axes. ........................ 35 Figure 2. 15: Ternary diagram to display focal mechanisms (Frohlich, 1992) ............. 36 Figure 2. 16: Difference between trend angle and plunge angle, N represents North. ............................................................................................................................ 37 Figure 2. 17: Stress and fault mechanisms .................................................................... 38 Figure 2. 18: Mohr diagram in 2-D with appropriate points in stress space ................. 41 Figure 2. 19: An example of Mohr circle within its associated stress shape ratio values (R) in an overall strike-slip regime. ................................................................... 45 Figure 2. 20: Different forces that give definition of friction. ...................................... 45 Figure 2. 21: Determination of the friction value from Mohr Coulomb failure criterion. ........................................................................................................................ 46 Figure 2. 22: Estimation of local magnitude using the Richter method. ...................... 48 Figure 2. 23: Illustration of the earthquake coda, which is normally used to estimate duration magnitude ......................................................................................... 50 Figure 2. 24: Frequency-domain transformation of the ground displacement wave. ... 52 Figure 2. 25: Deterministic seismic hazard analysis (Kramer, 1996) ........................... 57 Figure 2. 26: Probabilistic seismic hazard analysis. ..................................................... 58 Figure 2. 27: Source zone geometries (a) point source e.g. from a short fault, (b) linear source defined as shallow fault modeled, c) a three-dimensional source zone (Kramer, 1996). ............................................................................................................. 59 Figure 2. 28: Methods of source to site distance measurement. ................................... 59 Figure 2. 29: An example of curve showing the Gutenberg-Richter relation ............... 60 xxii Figure 2. 30: Kijko and Sellevoll (1992) method to obtain seismic hazard parameters.. ................................................................................................................... 66 Figure 2. 31: An example of logic tree in which different parameters are inserted for computing seismic hazard. ...................................................................................... 69 Figure 2. 32: Hypoellipse error ellipsoid relationships with SEH= MAXH/1.87, and SEZ=MAXZ/1.87. MAXZ is the radius of the 68% confidence ellipsoid (after Lahr, 1999). ................................................................................................................... 71 Figure 3. 1: Seismic stations in Madagascar ................................................................. 78 Figure 3. 2: Seismometers (1) CMG-3T, (2) STS-2, and (3) Trillium-120 .................. 81 Figure 3. 3: Schematic of a typical MACOMO seismic station. .................................. 82 Figure 3. 4: Installation of Vatomandry (VATO) station.. ........................................... 83 Figure 3. 5: A set of equipment inside the blue barrel (DAS – Power box- Batteries- Cables-PVC tube) ......................................................................................... 83 Figure 3. 6: Field servicing by Tsiriandrimanana Rakotondraibe and Fenitra Andriampenomanana in the SOLA station in northern Madagascar. ........................... 87 Figure 3. 7: Histogram showing data collected from MACOMO, local permanent, SELASOMA, and RHUM-RUM seismic stations. ....................................................... 89 Figure 3. 8: Characteristics of the waveform from signal-to-noise ratio that allowed reliable measurements of the peak amplitude and the spectral plateau for estimates of ML and MW ............................................................................................... 96 Figure 3. 9: Topographic map of Madagascar with 18,970 earthquakes from the combined MACOMO and Madagascar NDC catalogue between the time periods 1975-2016. .................................................................................................................... 97 Figure 3. 10: Earthquakes and seismic stations used to determine MW.. ...................... 97 Figure 4. 1: Seismicity map for Madagascar using data from the MACOMO project.. ........................................................................................................................ 100 Figure 4. 2: Computed depth distributions from (a) HYPOELLIPSE program (for 695 events) and (b) HYPODD program (for 422 events) ........................................... 101 Figure 4. 3: Earthquake distribution in Bemarivo area. .............................................. 103 xxiii Figure 4. 4: Earthquake distribution in Alaotra area. .................................................. 105 Figure 4. 5: Earthquake distribution in the Ankaratra area. ........................................ 106 Figure 4. 6: Earthquake distribution in the Itremo area (ITR). ................................... 107 Figure 4. 7: Map showing the earthquake distribution in Bekily area (BEK) from HYPOELLIPSE (grey circles) and the effects of event relocation by the HYPODD double-difference relocation method (red circles). ..................................................... 108 Figure 4. 8: Depth constraints for the 23 focal mechanisms, plotting the RMS travel-time residuals vs depth at intervals of 5 km. .................................................... 112 Figure 4. 9: Focal mechanisms of 23 earthquakes in Madagascar.. ............................ 114 Figure 4. 10: Focal mechanism solutions, constructed using P-wave first polarities and S-wave polarities, for 23 Madagascar intraplate earthquakes with good station coverage. ..................................................................................................................... 115 Figure 4. 11: Frequency-magnitude analysis of the 695 intraplate Madagascar earthquakes.. ................................................................................................................ 119 Figure 5. 1: A total of 12,727 earthquakes with ML>2.5 in Madagascar (2000- 2016) from the MACOMO bulletin and Madagascar NDC including all 182 earthquake events used for stress field inversion.. ...................................................... 127 Figure 5. 2: Focal mechanisms for 182 earthquake events ......................................... 128 Figure 5. 3: Ternary diagram showing the classification of the focal mechanisms used for the stress inversion calculation ..................................................................... 129 Figure 5. 4: Stress field inversion for 20 earthquakes in the Bemarivo zone (1). ...... 132 Figure 5. 5: Stress direction results for the Bemarivo area. ........................................ 133 Figure 5. 6: Same as Figure 5. 4 but for 23 events in Bekily zone. ....................... 134 Figure 5. 7: Stress direction results in the Bekily area. ............................................... 134 Figure 5. 8: Same as Figure 5. 4 but for the Famoizankova zone. ........................ 135 Figure 5. 9: Stress direction results in the Famoizankova region ............................... 136 Figure 5. 10: Same as Figure 5. 4 but for the events in the sub-zone A of Alaotra……………………. ........................................................................................ 137 xxiv Figure 5. 11: Same as Figure 5. 4 but for the events in the sub-zone B of Alaotra………………….. ........................................................................................... 138 Figure 5. 12: Same as Figure 5. 4 but for the events in the sub-zone C of Alaotra………………….. ........................................................................................... 139 Figure 5. 13: Stress direction results in the Alaotra region ......................................... 140 Figure 5. 14: Same as Figure 5. 4 but for the events in the sub-zone D of Ankaratra……………….. ........................................................................................... 142 Figure 5. 15: Same as Figure 5. 4 but for the events in the sub-zone E of Ankaratra……………….. ........................................................................................... 143 Figure 5. 16: Stress direction results in the Ankaratra region ..................................... 144 Figure 5. 17: Stress direction results in the Itremo region .......................................... 145 Figure 5. 18: Map showing the results of the minimum (blue) and maximum (red) stress direction after the inversion from each selected region. ................................... 147 Figure 6. 1: Plot of the estimated relation between Mw and ML for 80 eartthquakes reported in the MACOMO catalogue (black) and events derived from ISC bulletin (blue). .......................................................................................................................... 158 Figure 6. 2: Declustering using the Reasenberg (1985) method after testing different values into the parameter files. ..................................................................... 162 Figure 6. 3: Magnitude-time distribution and magnitude-distance distribution for declustered events in the whole area of Madagascar. ................................................. 164 Figure 6. 4: Timeline showing the annual number of earthquakes as a function of magnitude in the period 1975-2016, Mw≥2.9 .............................................................. 170 Figure 6. 5: Histogram showing the moment magnitude of all earthquakes in Madagascar, 1975-2016 inclusive. .............................................................................. 170 Figure 6. 6: Completeness estimation for all source zones for 42 years (1975- 2016, inclusive). Red stars indicate the MC values. .................................................... 171 Figure 6. 7: Results from ZMAP using data from the combined two catalogues (NDC+MACOMO) but plotted using GMT (Wessel and Smith, 1998). .................... 172 xxv Figure 6. 8: Madagascar’s location in the southeastern part of Africa https://platetectonicsafrica.weebly.com/cer.html ........................................................ 175 Figure 6. 9: Predicted ground motion for an M=3.5 earthquake with a 15 km focal depth calculated. .......................................................................................................... 176 Figure 6. 10: Map of Madagascar showing the 510 nodes (red squares) and seismic source zones (black boxes).......................................................................................... 178 Figure 6. 11: Logic tree for median ground motions for determining the branches at the final branches. ................................................................................................... 180 Figure 6. 12: Seismic hazard maps showing the expected PGA, estimated spectral acceleration for a period 1.0 s and 3.0 s at 10% probability of being exceeding in 50 year period (474-year return period). ..................................................................... 182 Figure 6. 13: Seismic hazard maps showing the expected PGA, estimated spectral acceleration for a period 1.0 s and 3.0 s at 2% probability of being exceeding in 50 year period (2,475-year return period). ....................................................................... 183 Figure 6. 14: Example of uniform hazard spectrum calculated for a site in the Ankaratra region in central Madagascar for a return period of 474- and 2,475- years (10% and 2% probability of exceedance in 50 years). ...................................... 184 Figure 6. 15: Earthquake hazard map for Africa in units of g with a 10% probability of being exceeded in 50 years (Poggi et al., 2017) ................................... 186 Figure 6. 16: Seismic hazard map showing PGA values in units of g with a 10% probability of being exceeded in 50 years (Razafindrakoto et al., 2009) ................... 187 xxvi List of tables Table 1. 1: Malagasy earthquakes with M≥5 ................................................................ 14 Table 1. 2: Crustal velocity models for Madagascar..................................................... 17 Table 1. 3: Focal mechanism solutions for Malagasy earthquakes. Source: (1) Shudofsky (1985), (2) Grimison and Chen (1988), (3) CMT-Harvard (http://www.globalcmt.org/CMTsearch.html ) and (4) Craig et al. (2011) .................. 18 Table 2. 1: Hypocentre quality factors based on the SEH and SEZ values (the vertical and horizontal errors) ....................................................................................... 71 Table 2. 2: Parameters for double difference relocation using hypoDD.inp ................ 73 Table 3. 1: Seismic stations used in the study (Shown in Figure 3. 1). ........................ 84 Table 3. 2: Stations that failed or were vandalized, 2011-2013. ................................... 87 Table 4. 1: Focal mechanism solutions for 23 Madagascar intraplate earthquakes .... 116 Table 4. 2: The numbers of earthquakes in the different binned ML magnitudes. ...... 120 Table 5. 1: Stress inversion results for 182 events recorded from 1 January 2000 - 31 December 2016 in Madagascar. ............................................................................. 146 AZ and PLG indicate Azimuth and Plunge. AZ is measured clockwise from the north and PLG downwards from the horizontal plane. ............................................... 146 Table 6. 1: Earthquakes used to determine the ML-Mw relation. (A) Earthquakes from the MACOMO bulletin, (B) Earthquakes from the ISC bulletin; “σ” is the average uncertainty value for Mw in m.u. (magnitude units). ..................................... 155 Table 6. 2: Parameters set for the declustering method following Reasenberg (1985) method after inserting 18,970 seismic events.................................................. 159 Table 6. 3: Parameter values inserted for declustering method using the Gardner and Knopoff (1974) approach. .................................................................................... 165 Table 6. 4: Mc values from the two bulletins NDC and MACOMO. (a) 1897-1974, (b) 1975-2008; (c) 2009-2013; (d) 2013-2016 ............................................................ 169 Table 6. 5: Largest earthquakes between November 1897 and December 1995 ........ 169 xxvii Table 6. 6: Characteristics of all seismic hazard parameters for all sources, where mmax is the maximum magnitude value, Beta (β) is equal to b*ln(10) a Gutenberg- Richter parameter. ....................................................................................................... 174 Table 6. 7: Parameters used to create the grid for seismic hazard .............................. 178 xxviii List of Acronyms BRSZ Bongolava-Ranotsara Shear Zone BRTT Boulder Real Time Technology CMAP Central Madagascar Alkaline Province CMT Centroid Moment Tensor DAS Data Acquisition System DSHA Deterministic Seismic Hazard Assessment FDSN Federation of Digital Seismograph Networks FPS Fault Plane Solution GFZ GeoForschungsZentrum (Geo-research Centre) GMT Generic Mapping Tools GPS Global Position System G-R Gutenberg-Richter HYPODD Hypo Double Difference IMAX Intensity Maximum IOGA ` Institut et Observatoire de Géophysique d'Antananarivo IRIS Incorporated Research Institutions for Seismology ISC International Seismological Center MACOMO Madagascar Comoros Mozambique MAG Magnitude Ma Million years from present Mya Million years ago MD Magnitude Duration ML Local Magnitude MW Moment magnitude NDC National Data Central NMAP Northern Madagascar Alkaline Province NRF National Research Foundation (South Africa) P Primary wave PASSCAL Portable Array Seismic Studies of the Continental Lithosphere PSHA Probabilistic Seismic Hazard Assessment xxix RHUM-RUM Réunion Hotspot and Upper Mantle – Réunions Unterer Mantel RMS Root Mean Square SMAP Southern Madagascar Alkaline Province SELASOMA Seismological Signatures in the Lithosphere/Asthenosphere system of Southern Madagascar SEH Standard Horizontal Errors SEZ Standard Vertical Errors 1 Chapter 1: General Introduction 1.1 Introduction The island of Madagascar was at the centre of a major rifting event that led to the break-up of Gondwana during the Jurassic and Cretaceous periods, which is indicated by its position relative to its one-time neighbors, East Africa and the Indian subcontinent. The seismicity of the island has been studied since the early decades of the 20th century (Poisson, 1924, 1930), while more detailed seismotectonic studies were initiated following the deployment of local seismic stations in 1977 (Bertil and Regnoult, 1998; Rindraharisaona et al., 2013). Clusters of earthquakes in the central part of Madagascar are believed to be related to volcanic activity during the Quaternary in the Ankaratra volcanic field, while linear zones of seismicity are attributed to active extensional faulting (Bertil and Regnoult, 1998; Norton and Sclater, 1979; Mahoney et al., 1991). The central region of Madagascar is the most seismically active. To date, no earthquake is known to have exceeded magnitude M6, and thus the level may be described as ‘moderate’. Recent deployments of seismic stations have helped to advance our understanding of the causes of seismic activity and the assessment of seismic hazard in Madagascar. In this study I analysed waveforms recorded by an island-wide temporary network of broadband stations, supplemented by waveforms from several permanent broadband and short-period stations. I used standard modern methods to locate hypocenters and determine fault plane solutions. I used the results from the central part of the island to define the orientation of the main stress axis. Finally, I carried out a Probabilistic Seismic Hazard Assessment (PSHA) in order to determine the spatial variation of Peak Ground Acceleration (PGA) values and the statistical value of earthquake occurrence in Madagascar. This preliminary chapter presents the problems to be solved and the principal objectives of this research following reviews of the geology of Madagascar and the results of previous geophysical studies. 2 1.2 Objectives This study uses high-quality newly available seismic data to examine the spatial distribution of earthquakes, fault characteristics, orientation of the stress field, and seismic hazard for the whole of Madagascar. Immediately prior to the MAdagascar-COmoros-MOzambique (MACOMO) project, 11 permanent stations were operational, four broadband and seven short-period. MACOMO deployed 28 temporary broadband seismic stations spread over the whole of Madagascar from 2011 to 2013 (details are given in Chapter 2). Standard methods were used to process, analyze and interpret that data (see Chapter 2). The novelty of this work lies in a significantly better understanding of the distribution and causes of seismicity in Madagascar, enabling a more reliable assessment of seismic hazard. In summary, this study aimed to: (a) Produce a new seismicity map for Madagascar using data from recently deployed temporary seismic stations, especially MACOMO; (b) Identify earthquake clusters and determine the focal mechanisms of selected earthquakes in order to define and characterize seismic source zones. (c) Determine the main direction of the stress field that is responsible for zones of active extensional faulting. (d) Assess the seismic hazard for Madagascar. 1.3 Outline This dissertation has three important topics, which are outlined in the following three sections. 1.3.1 Seismicity studies The first topic focuses on the distribution of earthquake activity during the period 2011-2013. Earthquakes recorded by the MACOMO network were located using the HYPOELLIPSE software (Lahr, 1989, 1999). The accuracy of the relative locations of nearby events was improved by applying the Double Difference 3 method implemented in HYPODD software (Waldhauser and Ellsworth, 2000; Waldhauser, 2001) (See Chapter 2). 1.3.2 Fault mechanisms and stress field model The second theme focuses on the seismotectonic structure of Madagascar. The geometry and form of the active faults in the central area of Madagascar was obtained after determining focal mechanisms using polarities of P- and S-wave first arrivals and the FOCMEC code (Snoke et al., 1984). 1.3.3 Magnitude scale and seismic hazard assessment For the final theme, an earthquake catalogue for the period 1975 to 2016 was compiled by merging data from bulletins produced by the National Data Centre (NDC), the International Seismological Centre (ISC) and the MACOMO project. Magnitude values were then homogenized to a common magnitude, MW, using a relation developed in this study. The seismic hazard parameters for identified source zones on the island were then determined and used in the seismic hazard calculation. 1.4 Problem statement Prior to this study, the seismotectonics and hence the seismic hazard in Madagascar was poorly characterized owing to shortcomings in the national seismic network. Most of the permanent seismic stations are single-component short-period instruments located in the central part of the island. There were some first-order questions that needed to be answered. For example: 1. Why are there volcanoes in the middle of this plate? 2. Does a boundary between the Nubian and Somalian plates (or various other micro-plates) pass through Madagascar? If so, where? 4 Even though the island-wide broadband MACOMO network was deployed for only two years, it gathered data of sufficient quantity and quality to attempt to answer these questions and produce an improved assessment of the seismic hazard. 1.5 Literature review: geological setting The current geological setting of Madagascar is due to two major rifting events: the separation of Madagascar from the Africa continent at ~160 Ma; and the separation of Seychelles and India from Madagascar at ~90 Ma (Kusky et al., 2007). The island once occupied a central position within the Gondwana supercontinent (Figure 1. 1, Stern, 1994), but rifted away from Africa (De Wit, 2003) and was displaced southwards parallel to the Davie Fracture Zone (Figure 1. 2 Coffin and Rabinowitz, 1992). However, its original position in Gondwana has been the subject of debate. Darracott (1974) and Förster (1975) place it in a position similar to that shown in Figure 1. 2 (i.e. it’s present day location); Wegener (1924) and Tarling (1971) placed Madagascar closer to present-day Mozambique; while others suggested that Madagascar was once adjacent to the Kenya-Somalia-Tanzania coast (e.g. Du Toit, 1937). Madagascar is divided into two major geological domains (Figure 1. 3). Precambrian rocks that pre-date Gondwana crop out in the eastern two-thirds of the island, while the western third is covered by younger sedimentary strata. Outcrops of younger volcanic rocks are found throughout Madagascar and range in age from the Cretaceous (beginning 145 Mya) to the Quaternary (see Figure 1. 3 and Figure 1. 4). 5 Figure 1. 1: Gondwana showing the position of Madagascar at the end of the Neoproterozoic Era (~550 Ma) (Kusky et al., 2003). Figure 1. 2: Displacement of Madagascar along the Davie Fracture Zone (Mascle et al., 1987). 6 1.5.1 Precambrian basement The Precambrian basement of Madagascar is bounded to the west by the Bongolava - Ranotsara shear zone (BRSZ) (Figure 1. 3). This Cambrian-age (542-488 Ma) structure was active during the last stages of the assembly of Gondwana, when the landmass consisting of present-day East Antarctica, East Africa and the Arabian-Nubian Shield coalesced, and is related to the Mozambican continental collision zone (Shackleton, 1996). The northern part of the Precambrian shield is composed of six geological units (Collins and Windley, 2002): Antongil Group, found in the northeastern region, inshore of Antongil Bay, where it forms the Antongil Craton; on the small islands of Nosy Boraha (Île Sainte Marie) and Nosy Mangabe; and in the central east coast, where it forms the Masora Craton ( i. Figure 1. 3). The group consists of migmatites, gneiss, and early Archaean granites, dated at about 3000 Ma (Jourde, 1971). The Antongil and Masora Cratons are Neoarchaean and Mesoarchaean in age with Indian (Dharwar Craton) affinities (Tucker et al., 1999; Kröner et al., 1999; Collins and Windley, 2002; Collins et al., 2012). ii. Tsaratanana Group, formed of gneiss, tonalites, ultramafic rocks, chromite-bearing paragneisses and greenstones, with zircon, dates between 2.75 and 2.49 Ga. (Tucker et al., 1999; Collins et al., 2000). The Tsaratanana and the Antananarivo blocks are divided by a mylonite zone (Tucker et al., 1999). iii. Bemarivo Group, defined by the Neoproterozoic amphibolite- and granulite-facies calc-alkaline volcanic, sedimentary and clastic rocks (Jöns et al., 2006; Thomas et al., 2008). The northern area of Madagascar is formed by the belt in the Bemarivo region. Low angle thrust faulting is observed in the southern part of the Bemarivo Belt. The Antongil Craton was thrust northwards over the Bemarivo block. iv. Antananarivo Group crops out in the central-northwest parts of the island. This Group consists of Neoarchaean gneisses, migmatites, granitoids and gabbros that were metamorphosed to granulite facies 7 between 700 and 550 Ma (Cox et al., 1998). Stratoid granites were emplaced at 630 Ma (Paquette and Nédélec, 1998). The east-west flexure that led to the formation of the north-south Angavo shear belt (Nédélec et al., 2003; Thomas et al., 2008) was accompanied by the intrusion of the porphyritic Carion granites at 532 Ma (Meert et al., 2001). v. Itremo Group is a Paleoproterozoic succession of low-grade shallow water sedimentary rocks composed of quartzite, mica schist and marble, and was subjected to low grade metamorphism. The Itremo Group was thrust onto the Antananarivo Group (Nédélec et al., 2003). vi. Ikalamavony Group, which is formed mostly of Proterozoic intermediate to mafic orthogneisses, formed between 1.03 Ga and 0.98 Ga, and has abundant volcano–clastic metasediments and quartzites (Tucker et al., 2011a, b). The southern part of the Precambrian shield is composed of three geological units (Hirtz et al., 1954; Boulanger and Pavlovsky, 1954; Bésairie, 1956, 1968-1971, Nicollet 1990). i. Vohibory Group consists of metasediments deposited during the Neoproterozoic Era from ~720 Ma to ~635 Ma (Collins et al., 2012) and amphibolite that are the products of granulite-facies metamorphism (Nicollet, 1990; Jöns and Schenk, 2008; Rakotovao, 2009). This Group also contains abundant mafic-ultramafic lenses of an unknown age. As the BRSZ has a left-lateral displacement, the Vohibory Group is tectonically associated with the eastern side of the Mozambique belt. ii. Ampanihy Group (also known as the Graphite Group) lies between the Androyen and Vohibory groups. It comprises a variety of rock types that have been metamorphosed to granulite facies. The most common rocks are graphite schists. 8 iii. Androyen Group consists of Neoproterozoic metasediments and graphitic gneiss. Figure 1. 3: Geological and tectonic map of Madagascar. Precambrian rocks and younger volcanics characterize the eastern and central parts of Madagascar, while sedimentary basins occupy the western third (after Bésairie, 1964). Dashed lines represent the shear zones: A - Ampanihy, AV - Angavo, BK - Bekily, BT - Betroka, BE - Betsileo, IF - Ifanadiana, T - Tranomaro, V - Vohibory. The red dashed line is the proposed plate boundary from Saria et al. (2014). 9 1.5.2 Sedimentary cover The sedimentary cover consists mostly of Karoo Supergroup rocks (Bésairie 1971, 1973; Hottin, 1972; See Figure 1. 3). The Palaeozoic-Cenozoic sedimentary deposits cover the western third of the island (the Antsiranana, Mahajanga and Morondava basins) and several very small areas close to the east coast (Figure 1. 3). In areas of Precambrian basement, sediments are limited to Cenozoic grabens such as the Lake Alaotra region and the Ankay graben, both probably Neogene in age (Bésairie and Collignon, 1971). 1.5.3 Volcanic fields Previous studies have shown that the island of Madagascar hosts two main volcanic fields (Besairie, 1973; Norton and Sclater, 1979) viz. Cretaceous-aged fields, which are mostly formed by tholeiitic basaltic lavas and Cenozoic-aged fields, formed by eruptions of alkali basalts. Cretaceous volcanic rocks are found along the periphery of the island. They are related to the Mesozoic rifting (ca. 95-85 Ma), when basaltic magmas enveloped all surface area of Madagascar (Storey et al., 1995). The Cretaceous volcanic rocks are preserved in the western sedimentary basin, along the eastern coast, and in the Volcan de l’Androy in the southern part of Madagascar (Figure 1. 3) (e.g. Mahoney et al., 1991; Rasamimanana, 1996). Three Cenozoic provinces are present (Figure 1. 3 and Figure 1. 4 ): i. The Northern Madagascar Alkaline Province (NMAP) consists of several volcanic fields. The Ambre-Bobaomby field, in the extreme northern part of the island, consists of Miocene volcanics capped by very recent cones. Important features include the Ankaizina volcano, cinder cones, lava flows, and several crater lakes (Lakes Anakarefo, Beandrazina and Mantsaboramadio, the Nanilezana Marsh and two other unnamed lakes). Two groups of cinder cones (viz. Nosy Be and 10 the Ankaizina field) are also found in the northern part of the island. These volcanoes are thought to have been active during the Holocene, although the date of the most recent eruption is not confirmed. Nosy Be, a volcanic island off the NW of Madagascar, contains < 1 Ma basaltic lava flows from well-preserved cinder cones. Several large crater lakes are observed in the central part of the island (Tucker and Moine, 2012). ii. The Central Madagascar Alkaline Province (CMAP) also consists of several fields. The Itasy field lies in central Madagascar (Figure 1. 4). The Kassigie volcano, situated about 100 km west of Antananarivo, is a prominent cone that covers an area of 100 km2. It lies directly on the Precambrian migmatite and gneiss, and formed from the Pliocene to the Holocene. The Ankaratra volcanic field covers a 5000 km2 area between Arivonimamo and Antsirabe (Figure 1. 3) and developed from the Miocene to the very recent Quaternary (Bésairie, 1973). iii. The Southern Madagascar Alkaline Province (SMAP). The most important field, the Ankililoaka, lies near the southwestern coast and is aged ca. 12 Ma (Cucciniello et al., 2018). 11 Figure 1. 4: Ankaratra and Lac Itasy volcanic fields in central Madagascar (after Bésairie, 1964). 1.6 Literature review: previous geophysical work Seismic activity in Madagascar is generally described as ‘moderate’ as few destructive earthquakes have been reported (Rindraharisaona et al., 2013). It was only by the end of the 1900s, when seismographs were installed on the island, that earthquake and structural seismology studies could be conducted, for example by Rakotondrainibe et al. (1977), Rambolamanana (1989), Bertil and Regnoult (1998), Razafindrakoto et al. (2009), and Rindraharisaona et al. (2013). Various topics have been investigated, such as the seismic velocity structure of the central area of Madagascar, focal mechanism studies, seismotectonic studies, and seismic hazard assessment. These studies were conducted using data from a network of 11 permanent seismic stations (Figure 1. 5). 12 1.6.1 Seismicity studies of Madagascar Historically, several large earthquakes have been recorded, notably an event with intensity VIII in the Itasy region in 1897, and with intensity VII in the Alaotra region in 1932. Earthquakes of magnitude M ≥ 5 recorded from 1897 to 1992 are listed in Table 1. 1. Bertil and Regnoult (1998) used data from the Institut et Observatoire de Géophysique d'Antananarivo (IOGA) to study the seismicity of Madagascar. They found that earthquakes in the central part of the island had an east-west extension focal mechanism and suggested that seismic activity is related to the volcanic activity. Rindraharisaona et al. (2013) found high levels of seismic activity in the Ankaratra, Alaotra and Ankay regions (Figure 1. 5). 13 Figure 1. 5: Seismicity and focal mechanisms determined by earlier studies.Red circles are events in the ISC Bulletin with magnitudes M>4 http://www.isc.ac.uk/iscbulletin/citing.php. Grey circles are events with magnitude M≥3 compiled by Rindraharisaona et al. (2013) using the Madagascar National Data Central (NDC) catalogue between 1990 and 2011. Fault plane solutions are from the Global Centroid Moment Tensor (CMT) catalogue (orange) (https://www.globalcmt.org/CMTsearch.html), and from previous work by Shudofsky (1985) (blue), Grimison and Chen (1988) (red), Foster and Jackson (1998) (grey), and from this study, MACOMO bulletin (black). http://www.isc.ac.uk/iscbulletin/citing.php https://www.globalcmt.org/CMTsearch.html 14 Table 1. 1: Malagasy earthquakes with M≥5 NDC catalogue (1897-1992) : ISC (1964-1992), intensities from Poisson (1930). DATE (dd/mm/yyyy) LAT LON MAG IMAX 03/11/1897 Itasy Region VIII 04/09/1932 Alaotra Region VII 29/03/1943 -18.6 44.2 ML 6.0 28/02/1945 Itasy Region VII 18/05/1965 -17.6 49.9 MB 5.4 30/06/1974 -13.47 48.90 MB 5.4 04/04/1975 -21.24 45.12 MB 5.3 03/01/1983 -23.50 44.60 MB 5.0 27/12/1983 -18.30 45.10 MB 5.1 04/10/1985 -18.30 48.40 MB 5.1 21/04/1991 -18.35 46.37 MB 5.8 29/01/1992 -13.56 49.93 MB 5.0 20/02/1992 -23.24 44.55 MB 5.0 14/11/1992 22.87 45.74 MB 5.1 1.6.2 Magnitude scale The earthquake magnitude scale for Madagascar was not well determined due to the shortage of earthquake data. The National Data Centre catalogue was compiled from a variety of sources. According to Rakotondrainibe (1977), the duration magnitude scale was developed using selected data recorded by the permanent short period stations and is in accord with the USGS mb for M<5.0. It is defined by: 𝑀𝐷 = 2.04 log 𝜏 + 0.0004 ∆ − 0.74 (1.1) where MD represents the magnitude duration, τ is the signal duration and 𝚫 the epicentral distance. 15 It was subsequently modified by Bertil and Regnoult (1998): 𝑀𝐷 = 2.04 log 𝜏 + 0.0004 ∆ − 0.74 + (1.75 x 10−7)𝜏2 (1.2) Magnitude duration MD and local magnitude ML were considered to be complementary for earthquakes having M<4.0, with a variance smaller than 0.2. 1.6.3 Velocity structure Several structural seismology studies have been carried out since the late 1900s. Rakotondrainibe (1977) used the origin and arrival times of earthquakes to compute, using a least-squares method, the apparent velocities of P- and S-waves in Madagascar (Table 1. 2, model 1). Rambolamanana et al. (1997) studied the crustal structure in the central area of Madagascar applying the simultaneous inversion of earthquake source parameters (Table 1. 2, model 2). Rindraharisaona et al. (2013) studied the lithospheric structure and seismicity of Madagascar. They estimated the characteristics of the structure underneath the central area of Madagascar using jointly inversion of Rayleigh wave group velocities and receiver functions (Table 1. 2 - model 3) and concluded that the lithosphere beneath Madagascar is thinner than that beneath East Africa. They used their new velocity model to locate earthquakes and determine focal mechanisms. We determined an ‘average’ 1D velocity model for Madagascar using results from Andriampenomanana et al. (2017) (Table 1. 2, model 4). The average 1D velocity model was derived from the VS-depth functions at each of the MACOMO stations reported by Andriampenomanana et al. (2017) (Figure 1. 6a) which were weighted according to the emplacement of the earthquake stations and the mapped geology (shown in Figure 1. 3) to produce an ‘average’ VS model for the whole island (See Appendix A). From the VP/VS ratios of 1.81 for the mantle and 1.75 for the crust, the corresponding VP model (Figure 1. 6) was determined. The validity of using a ratio of 1.75 was confirmed by a Wadati plot, while the validity of applying 1.81 for the mantle is confirmed by a Pn tomographic study (Andriampenomanana, 2017; Andriampenomanana et al., 2020). 16 a b c Figure 1. 6: (a) S-wave velocity model from all stations (Andriampenomanana et al. 2017), (b) Calculated vertical P-wave velocity structure for Madagascar using results from the (c) Wadati plot. 17 Table 1. 2: Crustal velocity models for Madagascar Depth (km) VP (km/s) VS (km/s) MODEL 1: Rakotondrainibe (1977) 0 5.9 3.4 10 5.9 3.4 12.5 6.7 3.8 25 6.7 3.8 30 6.7 3.8 35 8.1 4.6 40 8.1 4.6 45 8.1 4.6 50 8.3 4.7 MODEL 2: Rambolamanana et al. (1997) 0 5.9 3.4 10 6.5 3.7 20 6.7 3.8 30 6.3 3.6 35 6.3 3.6 42 8.1 4.6 50 8.1 4.6 MODEL 3: Rindraharisaona et al. (2013). 0 5.4 3.1 5 5.6 3.2 10 5.9 3.4 20 6.5 3.7 25 6.7 3.8 35 6.7 3.8 42 7.9 4.5 50 7.7 4.4 MODEL 4: Derived from Andriampenomanana et al. (2017) 0 5.6 3.2 10 5.9 3.4 13 6.4 3.6 16 6.5 3.7 30 7.0 4.0 33 7.2 4.1 40 7.9 4.4 55 7.9 4.4 18 1.6.4 Focal mechanisms Focal mechanisms could only be determined once sufficient permanent stations had been installed. The earliest earthquake for which a mechanism has been determined occurred in 1965 and was found to be extensional (Grimison and Chen, 1988). The predominant normal faulting mechanism was confirmed by five mechanisms determined by the global CMT study and one by Craig (2011). A large number of extensional mechanisms were identified in the central area of the island, which is in accordance with the earlier findings of Bertil and Regnoult (1998). The most recent study by Rindraharisaona et al. (2013) concetrated on focal mechanisms determination for Madagascar, she used earthquake data recorded by the permanent stations. A large number of extensional mechanisms were identified in the central area of the island, which is in accord with the earlier findings of Bertil and Regnoult (1998). Table 1. 3: Focal mechanism solutions for Malagasy earthquakes. Source: (1) Shudofsky (1985), (2) Grimison and Chen (1988), (3) CMT-Harvard (http://www.globalcmt.org/CMTsearch.html ) and (4) Craig et al. (2011) Source Date dd/mm/yyy Lat Lon Depth km Mw Strike Dip Rake 2 18/05/1965 -17.60 49.91 15 5.4 350 50 -95 2 04/04/1975 -21.24 45.13 13 5.6 95 75 44 1 04/04/1975 -21.24 45.13 11 5.4 245 85 140 3 31/01/1983 -23.00 44.39 15 5.1 326 45 -90 3 27/12/1983 -17.88 44.96 33 5.1 349 56 -84 3 04/101985 -18.08 48.62 10 5.5 357 61 -75 3 21/04/1991 -18.51 46.42 15 5.5 327 44 -128 3 14/11/1992 -23.01 45.54 15 5.0 350 45 -90 4 14/11/1992 -23.11 45.87 26 5.0 350 45 -90 http://www.globalcmt.org/CMTsearch.html 19 1.6.5 Hazard studies Seismic hazard studies began with the determination of the frequency-magnitude relationship. Earthquake data between the years 1979 to 1994 were used to establish an empirical relationship between magnitude (M) and the cumulative number of events greater than M for Madagascar (N) (Bertil and Regnoult, 1998). 𝑙𝑜𝑔𝑁 = 9.23 − (1.65 ± 0.10)𝑀 (1.3) The b-value of 1.65 value is high; values close to 1 are generally produced by natural tectonic seismicity (as opposed to induced seismicity, where values may range from 0.5 to 1.5). This may be due to the small number of earthquakes in the catalogue and also the fact that only short period permanent stations were used to calculate the magnitudes. This value could imply uncertainties of the interpretation results. Razafindrakoto et al. (2009), using a neodeterministic approach which followed the calculation of realistic synthetic seismograms to estimate regional ground motion, calculated peak acceleration, velocity and displacement values of 0.03 g, 1.6 cm/s and 0.5 cm, respectively, for a site in Madagascar. They concluded that these values were below the damage threshold and consequently that the risk posed by earthquakes was low. 20 Chapter 2: Theory In this chapter I describe the theory that underlies the techniques that I applied to achieve the main objectives stated in Chapter 1. 2.1 Seismic waves Seismology is the study of earthquakes (and other sources of ground vibration), the propagation of the elastic waves that they generate and the shaking that they produce. Earthquakes may be caused by natural phenomena such as fault rupture and volcanic activity, as well as by human activities such as mining and the impoundment of large reservoirs. Seismograms may also be analysed to deduce the structure of the Earth’s interior. Seismic waves produced by earthquakes are grouped into two principal categories: surface waves that are guided by the Earth’s surface and body waves that travel through Earth’s interior (Figure 2. 1).There are two types of body waves, viz. compressional waves (P-wave or primary wave), which consist of alternating expansion and compression pulses in the direction of the wave displacement; and shear waves (S-wave or secondary wave), whose displacement is perpendicular to the direction of movement of the wave (Figure 2. 1). 21 Figure 2. 1: Main types of waves produced by an earthquake. (https://courses.lumenlearning.com/geo/chapter/reading-studying-the- earths-interior ) Figure 2. 2 shows a seismogram recorded by the local permanent seismic station SBV of a magnitude 3.5 earthquake on 16 March 2012 to illustrate typical P and S waves as they were recorded. Figure 2. 2: Vertical component seismogram recorded at the SBV seismic station produced by a M3.5 earthquake on March 16th, 2012 at 11:56:55 UTC. P and S mark the onsets of the primary or P-wave and secondary or S-wave, respectively. https://courses.lumenlearning.com/geo/chapter/reading-studying-the-earths-interior https://courses.lumenlearning.com/geo/chapter/reading-studying-the-earths-interior 22 2.2 Source parameters and mechanisms 2.2.1 Location and origin time The focus of an earthquake (also called the hypocenter) is the point within the Earth where the rupture initiates, while the epicenter is the point on the Earth’s surface directly above it (Figure 2. 3). An earthquake is generally defined as ‘shallow’ if the focus is at a depth of <100 km from the ground surface, ‘intermediate’ if the focus is at depths between 100 km and 300 km, and ‘deep’ if the focus is >300 km. The greatest depth at which earthquakes occur is about 700 km. Figure 2. 3: Earthquake focus (hypocenter) and epicenter. The ground surface and subsurface are colored green and grey, respectively. (https://earthquakesandplates.files.wordpress.com/2008/05/epicenter.gif ) The first task in seismogram analysis is the identification and picking of the first phase arrivals, then measurement of the amplitude of ground motion, followed by the determination of the time, epicenter and depth of the point where the earthquake rupture initiated. In principle, the earthquake location can be estimated using three or more seismic stations using the ‘triangulation method’ (Figure 2. 4). Firstly, the distance from each station to the focus is calculated using the difference between the arrival times of the P- and S-waves and prior knowledge of the average P- and S-wave velocities. Secondly, circles with radii equal to the focal distance are drawn around each station. Lastly, the location is determined https://earthquakesandplates.files.wordpress.com/2008/05/epicenter.gif 23 from the intersection of these circles. The solutions are more accurate if the Earth’s velocity structure is well-known and the refraction of waves is taken into account, and the stations are well-distributed. Figure 2. 4: Triangulation method to locate earthquakes. http://www.bbc.co.uk/education/guides/zq7ycdm/revision/2 ) Several computer programs using different mathematical algorithms have been elaborated to determine earthquake locations. In this study SEISAN (Havskov and Ottemöller, 2010), HYPOELLIPSE (Lahr 1989, 1999), and HYPODD (Waldhauser, 2000, 2001) were used. Many location methods use the iterative least square technique (Chiu, 2006), including the three codes that I used. An initial guess of the coordinates of the focus is made, usually the coordinates of the station that first recorded it and hence closest to it. The theoretical arrival times produced by an earthquake at this hypothetical location are then compared with the actual arrival times at each station. The difference is known as the travel time residual. The discrepancies between the theoretical and actual times are then used to adjust the location of the hypothetical focus. This process is repeated until there is no further reduction in the sum of the squares of the residuals. http://www.bbc.co.uk/education/guides/zq7ycdm/revision/2 24 Here we review the method developed by Geiger (1910, 1912). The coordinate of the ith station is (xi, yi, 0), representing that the station is on the surface, and the corresponding arrival time is tj. Let fi(x) be the arrival time function at the ith station, and x represents the hypocentral parameters and T the transpose of the column vector: x= (x, y, z, t)T The value of fj(xo) at a nearby location, xo, can be expressed by a first-degree Taylor polynomial: 𝑓𝑗 (𝒙) = 𝑓𝑗(𝒙0 + 𝛿𝒙) = 𝑓𝑗(𝒙0) + 𝜕𝑓𝑗 𝜕𝑥 𝛿𝑥 + 𝜕𝑓𝑗 𝜕𝑦 𝛿𝑦 + 𝜕𝑓𝑗 𝜕𝑧 𝛿𝑧 + 𝜕𝑓𝑗 𝜕𝑡 𝛿𝑡 (2.1) with x= x0+𝛿𝑥 (2.2) and x0 = (x0, y0, z0, t0) T (2.3) and 𝛿𝑥 = (𝛿𝑥, 𝛿𝑦, 𝛿𝑧, 𝛿𝑡)𝑇 (2.4) Equation (2.1) can be expressed in vector notation: 𝑓𝑗(𝑥0 + 𝛿𝑥) = 𝑓𝑗(𝑥0) + 𝑔𝑗 𝑇𝛿𝑥 (2.5) where 𝑔𝑗 𝑇 is the transpose of the gradient vector gj and defined as follows: 𝑔𝑗 𝑇 = 𝛻𝑓𝐽(𝒙) = ( 𝜕𝑓𝑗 𝜕𝑥 , 𝜕𝑓𝑗 𝜕𝑦 , 𝜕𝑓𝑗 𝜕𝑧 , 𝜕𝑓𝑗 𝜕𝑡 ) (2.6) 𝑓𝑗(𝒙0 + 𝛿𝒙) represents the arrival time recorded at the jth station (the observed arrival time) and 25 𝜕fj 𝜕𝑥 𝛿𝑥 + 𝜕fj 𝜕𝑦 𝛿𝑦 + 𝜕fj 𝜕𝑧 𝛿𝑧 + 𝜕fj 𝜕𝑡 𝛿𝑡 (2.7) Equation (2.7) is the correction value, which the partial derivatives function of the hypocentral parameters. For calculating a system of the Equation (2.5), it is important to find x0 such that the estimated arrival times will best fit the observed arrival times hence x0 can be defined as the hypocenter of the event. Equation (2.7) can be writen into the form: 𝜕𝑓𝑗 𝜕𝑋 𝛿𝑥 + 𝜕𝑓𝑗 𝜕𝑌 𝛿𝑦 + 𝜕𝑓𝑗 𝜕𝑍 𝛿𝑧 + 𝜕𝑓𝑗 𝜕𝑇 𝛿𝑡 = 𝑟𝑗 (2.8) where 𝑟𝑗 = 𝑡𝑜𝑏𝑠 − 𝑡𝑐𝑎𝑙𝑐 (2.9) 𝑡𝑜𝑏𝑠 = fj(x), 𝑡𝑐𝑎𝑙𝑐= fj(x0), and rj is the residual In matrix notation, Equation (2.8) can be written: 𝐴𝛿𝒙 = 𝒓 (2.10) where; A =[ 𝜕𝑓1 𝜕𝑋 𝜕𝑓1 𝜕𝑌 𝜕𝑓1 𝜕𝑍 𝜕𝑓1 𝜕𝑍 ⋮ 𝜕𝑓𝑚 𝜕𝑋 𝜕𝑓𝑚 𝜕𝑌 𝜕𝑓𝑚 𝜕𝑍 𝜕𝑓𝑚 𝜕𝑍 ]; (2.11) 𝛿𝒙 =[ 𝛿𝑥 𝛿𝑦 𝛿𝑧 𝛿𝑡 ] and r= [ 𝑟1 ⋮ 𝑟𝑚 ] The least squares solution to the system by Equation (2.10) satisfies (Strang, 1976) 26 𝐴𝑇𝐴𝛿𝒙 = 𝐴𝑇𝒓 (2.12) or 𝛿𝒙 = (𝐴𝑇𝐴)−1𝐴𝑇𝐫 (2.13) The total effect of the m mismatches between the observed and calculated arrival times is the event residual or simply residual; it is defined by the least-squares solution (Ge, 1995). 𝑅𝑒𝑠 = √ 𝑟𝑇 𝑟 𝑚 − 𝑞 (2.14) where m is the number of equations and q is the degree of freedom. The number of hypocentral parameters defined by the Equation (2.13) has degree of freedom 4. The correction vector, 𝛿𝒙, is determined and it can be introduced to the previous trial solution to form a new one. The iteration is carried out until the given error criterion was performed and the last trial solution is then considered as the solution for the hypocenter. 2.2.2 Correction for instrument response Mathematically, a seismogram is the combined product of the signal produced by the earthquake with impulsive response from the equipment that recorded the signals (Havskov and Ottemōller, 2010). The signal recorded from a seismic station is represented by the equation: 𝑥(𝑛) = ℎ(𝑛) ∗ [𝑠(𝑛) + 𝑧1(𝑛) + 𝑧2(𝑛)] (2.15) where x(n) represents the output of the seismograph, h(n) is the impulse response of the seismograph, s(n) is the actual ground motion produced by the earthquake, z1(n) is the natural noise, and z2(n) is the noise obtained from the acquisition system. After combining the different noises 𝑧(𝑛) = [ℎ(𝑛) ∗ 𝑧1(𝑛)] + [ℎ(𝑛) ∗ 𝑧2(𝑛)] Equation (2.15) can be written as: 𝑥(𝑛) = ℎ(𝑛) ∗ 𝑠(𝑛) + 𝑧(𝑛) (2.16) 27 In the frequency domain, multiplication is the equal to the convolution operation. Thus Equation (2.16) can be written as: 𝑋(𝑓) = 𝐻(𝑓). 𝑆(𝑓) + 𝑍(𝑓) (2.17) where X(f), H(f), S(f) and Z(f) are the Fourier Transforms of the seismogram, impulse response, actual ground motion and noise, respectively. It is therefore necessary to apply a deconvolution method to remove the effect of the instrument from the seismogram. Deconvolution can be carried out in either the time or the frequency domains. In the frequency domain 𝑆(𝑓) = 𝑋(𝑓) 𝐻(𝑓) − 𝑍(𝑓) 𝐻(𝑓) (2.18) 2.2.3 Fault geometry When an earthquake occurs, the amplitude and polarity of the radiated seismic waves are a product of the fault geometry and the slip direction. The fault geometry is characterized by the strike angle (𝜙𝑓), dip angle (δ), and slip angle (λ). The strike direction is the intersection of the fault plane with the Earth’s surface and varies from 0° to 360°, measured clockwise from North; the dip angle is the angle between the horizontal plane and fault plane and varies from 0° to 90°; and the slip angle varying from +180° to -180° from the horizontal on the fault plane. d is the slip direction, n is the normal to the fault plane, and x1, x2 and x3 are the Cartesian coordinate system axes (Figure 2. 5; Figure 2. 6). The strike is characterized by the equation: tan[180° − (𝜙𝑓 + 90°)] = − 𝑥2 𝑥1 (2.19) where x1, x2 are defined in Figure 2. 14 and 𝜙𝑓 = arctan (− 𝑥1 𝑥2 ) (2.20) Dip angle is defined by: 28 sin(90° − 𝛿) = sin(𝛿 − 90°) = 𝑐𝑜𝑠𝛿 = − 𝑥3 |𝑛| (2.21) with δ= arccos(-x3) Rake is defined by: sin(180° − 𝜆) = − sin(𝜆 − 180°) = −𝛾3/𝑠𝑖𝑛𝛿 |𝑑| (2.22) where 𝜸𝟑 𝐚𝐧𝐝 𝜹 are shown in Figure 2. 6, and 𝜆 = arcsin(− 𝛾3 𝑠𝑖𝑛𝛿 ) (2.23) and ( −sin 𝛿 sin𝜙𝑓 𝑠𝑖𝑛𝛿𝑐𝑜𝑠𝜙𝑓 −𝑐𝑜𝑠𝛿 ) = 𝒏 (2.24) ( 𝑐𝑜𝑠 𝜆𝑐𝑜𝑠𝜙𝑓 + cos 𝛿𝑠𝑖𝑛𝜆𝑠𝑖𝑛𝜙𝑓 cos𝜆 𝑠𝑖𝑛𝜙𝑓 − cos𝛿 sin 𝑐𝑜𝑠𝜙𝑓 −sin𝛿 𝑠𝑖𝑛 ) = 𝜸 (2.25) Figure 2. 5: Geometrical description of a fault plane 𝜸 is the rake direction and n is the normal to the fault plane (http://www.geol.tsukuba.ac.jp/~yagi- y/text/Source_meca_v0.7.pdf). http://www.geol.tsukuba.ac.jp/~yagi-y/text/Source_meca_v0.7.pdf http://www.geol.tsukuba.ac.jp/~yagi-y/text/Source_meca_v0.7.pdf 29 a) Strike b)Dip c)Rake Figure 2. 6: Illustration of the method for calculating strike, dip, and rake angles (http://www.geol.tsukuba.ac.jp/~yagi-y/text/Source_meca_v0.7.pdf) 2.2.4 Focal mechanism solutions Focal mechanism solutions use the ground motions produced by an earthquake and observed at many seismic stations to estimate the geometry of the fault plane (i.e. strike and dip) and the direction of the slip vector. Methods used to determine focal mechanism solutions include P-wave first motion polarity (Stein and Wysession, 2003), P/S amplitude ratios analysis (Kisslinger, 1980), polarization and amplitude of S-waves (Khattri, 1973), and moment tensor inversion. In this study we determined focal mechanism solutions using the polarity of the first-motions of the P- and S-waves recorded at each station. The polarity (compressive or dilatational in the case of P-waves; or direction of the first motion of the SV or SH pulse, in the case of S-waves) is plotted on a stereogram at the point where the ray emanating from the focus pierces the focal sphere (a small imaginary sphere centered on the focus). The observed polarities are compared with those produced by a systematic search of all possible fault plane and slip vector orientations, and a set of plausible solutions is determined. The position on the focal sphere where the polarity of ground motion is plotted depends on the azimuth (Φ) and take-off angle. The azimuth is simply the angle http://www.geol.tsukuba.ac.jp/~yagi-y/text/Source_meca_v0.7.pdf 30 describing the direction of the line joining the epicenter and station, measured clockwise from the North direction Figure 2. 7. Figure 2. 7: Measurement of the azimuth angle https://service.iris.edu/irisws/rotation/docs/1/help/ The take-off angle is the angle at which the seismic ray leaves the earthquake source during the signal’s transference from the earthquake focus to the station through the Earth. It is more complicated to determine. Consider, for instance, two stations lying on the same azimuth but at different epicentral distances. The ray path to each station will differ, depending on the focal depth, the Earth’s velocity structure, and the distance to the station. A standard Earth model is used to calculate the take-off angle (Figure 2. 8). It can also be determined from the gradient of the travel-time curve. In general, the greater the epicentral distance, the smaller the take-off angle. Figure 2. 8: Definition of the take-off angle https://nctr.pmel.noaa.gov/education/ITTI/seismic/Focal_mech_USGS.pdf https://service.iris.edu/irisws/rotation/docs/1/help/ https://nctr.pmel.noaa.gov/education/ITTI/seismic/Focal_mech_USGS.pdf 31 In the P-wave polarity method, focal mechanisms are determined following three steps: (i) determine the location of the earthquake (latitude, longitude, depth), (ii) Determine the azimuth and take-off angle values, and (iii) Determine and plot first-arrival P-wave polarities (an upward inflection on the vertical component indicates a compression, while a downward inflection indicates a dilatation; (see Figure 2. 9), plot on the stereonet, and find the nodal planes that best separate the observations into quadrants. Figure 2. 9: Illustration of P-wave upward and downward first motions (Havskov and Ottemöller, 2010) The picking of S-wave polarities is more complicated: 1. The two horizontal components (E-W and N-S) are rotated into radial (R) and transversal (T) components (Figure 2. 10 a). 2. A second rotation into the 3D LQT ray system is then performed Figure 2. 10 b). The R and T components are combined with the vertical (Z) component to determine the motion in the longitudinal (L) direction (along the ray path) and the Q-direction (perpendicular to the ray path in the vertical plane). 3. The polarities of the ground motion in the Q-direction (SV) and the T- direction (SH) are then determined. 32 (a) (b) Figure 2. 10: Transformation of ZNE to ZRT and LQT (https://service.iris.edu/irisws/rotation/docs/1/help/) Transformation from ZNE to ZRT is as follows: [ 𝑅 𝑇 𝑍 ] = [ cos 𝜃 sin 𝜃 0 − sin 𝜃 cos 𝜃 0 0 0 1 ] [ 𝐸 𝑁 𝑍 ] with 𝜃 = 3𝜋 2 − 𝜉 (2.26) 𝜉 is the back azimuth that represents the vector pointing from the seismic station to the epicenter, measured clockwise from the north (Figure 2. 10). The transformation from ZNE to LQT is given by: [ 𝐿 𝑄 𝑇 ] = [ cos 𝜃 sin 𝛿 sin 𝜉 − sin𝛿 𝑐𝑜𝑠 𝜉 sin 𝜃 cos𝛿 sin 𝜉 cos𝛿 cos 𝜉 0 −𝑐𝑜𝑠𝜉 𝑠𝑖𝑛𝜉 ] [ 𝑍 𝐸 𝑁 ] (2.27) with 𝛿 the angle of the incident wave from the vertical, 𝛿 = 0𝑜 for vertically incident wave 𝛿 = 90𝑜 for horizontally incident wave Example of picks from P-wave and S-wave are presented in Figure 2. 11. https://service.iris.edu/irisws/rotation/docs/1/help/ 33 Figure 2. 11: Examples of polarity picking on a local M3.1 event recorded at station BKTA on 4 July 2013 at 13:40:59.08. P-wave picks on the (a) E- W, (b) (N-S) and (c) Z component; S-wave picks after rotating from ZNE to ZRT on the (d) R, (e) T and (f) Z component. The advantage of picking both P- and S-wave polarities is to constrain the characteristics of the fault solutions better. In this work, we tried to pick both P- and S-waves polarities, although, due to the waveform characteristics some S- wave polarities were difficult to identify. In general, faulting can be grouped into three types: strike-slip (left lateral or right lateral), normal and reverse (See Figure 2. 12). The P-wave first motions are divided by two nodal planes into four quadrants: two opposite quadrants that experience compressional first motions, and two opposite quadrants that experience dilatational first motion. One nodal plane is the fault plane, while the other is termed the auxiliary plane (Figure 2. 13). Other observations are needed to identify which nodal plane is the fault plane e.g., surface rupture, aftershock distribution. By convention, the compressive quadrants are colored black (Figure 2. 14). The exact appearance of the stereogram will depend on the orientation of the planes and the slip vector. These diagrams are often referred to as ‘beach balls’. (a) (f) (c) (b) (a) (e) (d) 34 a) Normal faulting b) Reverse faulting b) Left-lateral strike-slip c) Right-lateral strike slip Figure 2. 12: Fault mechanisms and ‘beach ball’ diagrams. Black and white represent the compressional and dilatational quadrants, respectively. https://inventionsky.com/normal-reverse-and-strike-slip-faults/ https://inventionsky.com/normal-reverse-and-strike-slip-faults/ 35 Figure 2. 13: Focal mechanism for a normal faulting earthquake with the actual fault plane and auxiliary plane in a section BB'. The two nodal planes are the auxiliary plane and the fault plane (http://www- udc.ig.utexas.edu/external/TXEQ/faq_basics.html) From the focal mechanism, the orientations of the P-and T-axes that relate to the white and dark segments of the beach ball are determined (Figure 2. 14). P- and T- axes are geometric axes, not directly stress axes. They represent the average orientation of the quadrants in which the extensional and compressional stress axes lie. The P-axis, orthogonal to the B- and T-axes, is the pressure or compressional axis that bisects the dihedral angle between the two nodal planes. The T-axis is orthogonal to both B- and P-axes. The B-axis is the neutral axis, which is the intersection of the two nodal planes and is orthogonal to the P- and T- axes. Figure 2. 14: Normal fault plane solution showing the P- and T-axes https://www.usgs.gov/natural-hazards/earthquake-hazards/science/focal- mechanisms?qt-science_center_objects=0#qt-science_center_objects. Earthquake focal mechanisms can be displayed in a triangular diagram named the ‘ternary diagram’ (Frohlich, 1992; Frohlich and Apperson, 1992). The triangle http://www-udc.ig.utexas.edu/external/TXEQ/faq_basics.html http://www-udc.ig.utexas.edu/external/TXEQ/faq_basics.html https://www.usgs.gov/natural-hazards/earthquake-hazards/science/focal-mechanisms?qt-science_center_objects=0#qt-science_center_objects https://www.usgs.gov/natural-hazards/earthquake-hazards/science/focal-mechanisms?qt-science_center_objects=0#qt-science_center_objects 36 classifies the mechanisms as normal, thrust or strike-slip in terms of the dip angles or plunge of the P-, T-, and B-axes (Figure 2. 15). Figure 2. 15: Ternary diagram to display focal mechanisms (Frohlich, 1992) The diagram displays the focal mechanism in terms of the plunge angles for the P- , T- and B-axes (Figure 2. 16). Angles are notated as dip-P (δP), dip-T (δT), and dip-B (δB.). These values satisfy the equation 𝑠𝑖𝑛2 𝛿𝑇 + 𝑠𝑖𝑛2 𝛿𝑃 + 𝑠𝑖𝑛2 𝛿𝐵 = 1 (2.28) By considering the vertical of the three axes, 𝑥 = sin 𝛿𝑇 𝑦 = sin 𝛿𝑃 𝑧 = sin 𝛿𝐵 (2.29) and the Equation 2.28 becomes; 𝑥2 + 𝑦2 + 𝑧2 = 1 (2.30) The horizontal hG and the vertical distances vG from the center of the triangle (G) are defined by: ℎ𝐺 = 𝑐𝑜𝑠𝛿𝐵𝑠𝑖𝑛 𝜓/[(𝑠𝑖𝑛𝜃0. sin 𝛿𝐵 + 𝑐𝑜𝑠𝜃0𝑐𝑜𝑠𝛿𝐵𝑐𝑜𝑠 𝜓)] 𝑣𝐺 = (𝑐𝑜𝑠𝜃0𝑠𝑖𝑛𝛿𝐵 − 𝑠𝑖𝑛𝜃0𝑐𝑜𝑠𝛿𝐵 cosψ)/[(𝑠𝑖𝑛𝜃0𝑠𝑖𝑛𝛿𝐵 + 𝑐𝑜𝑠𝜃0𝑐𝑜𝑠𝛿𝐵𝑐𝑜𝑠𝜓)] (2.31) where 37 𝜓 is defined by 𝛹 = 𝑡𝑎𝑛−1 ( 𝑠𝑖𝑛𝛿𝑇 𝑠𝑖𝑛𝛿𝑃 ) − 45𝑜 𝜃0 = 35.26° is the dip angle of the P- T- and B-axes for the source mechanism that plots in the exact centre of the ternary diagram Figure 2. 16: Difference between trend angle and plunge angle, N represents North. In this study, the ternary diagram plot was implemented in MATLAB using the Frohlich (1992) algorithm. 2.3 Stress field 2.3.1 Basic concepts Stress describes the force acting on body and is expressed as the force per unit area. The stress vector may be decomposed into the component acting normal to a plane, which may be either compressive or tensile; and the component acting parallel to a plane, known as the shear stress. The stress state within a body is often described in terms of a set of orthogonal ‘principal stresses’ that act on orthogonal planes such that they do not produce any shear stress on them. They are designated 1, 2 and 3, indicating the maximum, intermediate and minimum stress, respectively. The stress field may also be described in terms of the vertical, minimum, and maximum horizontal stresses V, Hmin, and Hmax respectively. An earthquake may occur if the shear stress acting on a geological weakness exceeds its shear strength (Figure 2. 17). High horizontal stresses (i.e. 1 sub- 38 horizontal; Hmax>Hmin>V) will generally cause folding, reverse or thrust faulting, and the crust may thicken. If the horizontal stresses are relatively low (i.e. 1 sub-vertical; V>Hmax>Hmin), the dominant mode of faulting will be normal and the crust may thin. The third case (i.e. 1 sub-horizontal; Hmax>V>Hmin) will likely give rise to strike-slip faulting. There are, of course, many intermediate cases, which give rise to oblique-slip. Knowledge of the stress field is useful when analyzing the distribution of seismicity, the kinematics of continental deformation and assessing seismic hazard. The stress field may be analysed using knowledge of fault plane geometries, Mohr circle diagram, and tensor calculus. Figure 2. 17: Stress and fault mechanisms. http://www.geosci.usyd.edu.au/users/prey/Teaching/Geol- 3101/Strain/stress.html 2.3.2 Mohr’s circle of stress The Mohr circle diagram gives the relationship between the orientation of any plane and the values of the normal and shear stress acting on it (Holtz et al., 1981). This diagram can show the stress in two and three dimensions. In this study, Mohr’s circle diagram is presented into three dimensions. According to Cauchy’s http://www.geosci.usyd.edu.au/users/prey/Teaching/Geol-3101/Strain/stress.html http://www.geosci.usyd.edu.au/users/prey/Teaching/Geol-3101/Strain/stress.html 39 law and with respect to the principal direction 1, 2, and 3; the normal stress is defined by: 𝜎𝑁 = 𝑡 (𝑛). 𝑛 = (𝜎𝑛). 𝑛 (2.32) where 𝑡(𝑛) is the traction vector, 𝜎𝑛 represents the normal stress on the plane with normal n. In term of matrix calculus, this can be written: 𝜎𝑁 = {[ 𝜎1 0 0 0 𝜎2 0 0 0 𝜎3 ] [ 𝑛1 𝑛2 𝑛3 ]} [ 𝑛1 𝑛2 𝑛3 ] = 𝜎1 𝑛1 2 + 𝜎2 𝑛2 2+𝜎3 𝑛3 2 (2.33) Since 𝑛1 2 + 𝑛2 2 + 𝑛3 2 = 1 with 𝜎1 ≥ 𝜎2 ≥ 𝜎3 𝜎1 = 𝜎1 (𝑛1 2 + 𝑛2 2 + 𝑛3 2) ≥ 𝜎1 𝑛1 2 + 𝜎2 𝑛2 2+𝜎3 𝑛3 2 = 𝜎𝑁 (2.34) Similarly, 𝜎𝑁 = 𝜎1 𝑛1 2 + 𝜎2 𝑛2 2+𝜎3 𝑛3 2 ≥ 𝜎3 (𝑛1 2 + 𝑛2 2 + 𝑛3 2) ≥ 𝜎3 (2.35) As a result, the maximum normal stress occurring in a point is the maximum principal stress and the minimum normal stress is the minimum principal stress. The shear stress occurring on the plane is given by: 𝜎𝑠 2 = (𝜎1 2𝑛1 2 + 𝜎2 2𝑛2 2 + 𝜎3 2𝑛3 2) − (𝜎1 𝑛1 2 + 𝜎2 𝑛2 2+𝜎3 𝑛3 2 )2 (2.36) As 𝑛1 2 + 𝑛2 2 + 𝑛3 2 = 1, we can write 𝜎𝑠 2 + 𝜎𝑁 2 = 𝜎1 2𝑛1 2 + 𝜎2 2𝑛2 2 + 𝜎3 2𝑛3 2 (2.37) So we have the system of equations: 𝜎𝑠 2 = (𝜎1 2𝑛1 2 + 𝜎2 2𝑛2 2 + 𝜎3 2𝑛3 2) − (𝜎1 𝑛1 2 + 𝜎2 𝑛2 2+𝜎3 𝑛3 2 )2 𝑛1 2 + 𝑛2 2 + 𝑛3 2 = 1 𝜎𝑠 2 + 𝜎𝑁 2 = 𝜎1 2𝑛1 2 + 𝜎2 2𝑛2 2 + 𝜎3 2𝑛3 2 (2.38) Solving these equations produces: 𝑛1 2 = (𝜎𝑁 − 𝜎2 )(𝜎𝑁 − 𝜎3 ) + 𝜎𝑠 2 (𝜎1 − 𝜎2 )(𝜎1 − 𝜎3 ) 𝑛2 2 = (𝜎𝑁 − 𝜎3 )(𝜎𝑁 − 𝜎1) + 𝜎𝑠 2 (𝜎2 − 𝜎3 )(𝜎2 − 𝜎1 ) (2.39) 40 𝑛3 2 = (𝜎𝑁 − 𝜎1)(𝜎𝑁 − 𝜎2 ) + 𝜎𝑠 2 (𝜎3 − 𝜎1 )(𝜎3 − 𝜎2 ) By taking 𝜎1 ≥ 𝜎2 ≥ 𝜎3with the condition that the squares of the normal components must be positive, one has that (𝜎𝑁 − 𝜎2 )(𝜎𝑁 − 𝜎3 ) + 𝜎𝑠 2 ≥ 0 (𝜎𝑁 − 𝜎3 )(𝜎𝑁 − 𝜎1) + 𝜎𝑠 2 ≤ 0 (𝜎𝑁 − 𝜎1)(𝜎𝑁 − 𝜎2 ) + 𝜎𝑠 2 ≥ 0 (2.40) and we can get 𝜎𝑠 2 + [𝜎𝑁 − 1 2 (𝜎2 + 𝜎3)] 2 ≥ [ 1 2 (𝜎2 − 𝜎3)] 2 𝜎𝑠 2 + [𝜎𝑁 − 1 2 (𝜎1 + 𝜎3)] 2 ≥ [ 1 2 (𝜎1 − 𝜎3)] 2 𝜎𝑠 2 + [𝜎𝑁 − 1 2 (𝜎1 + 𝜎2)] 2 ≥ [ 1 2 (𝜎1 − 𝜎2)] 2 (2.41) Considering the coordinates (𝜎𝑁,𝜎𝑆), a plot of circles in (𝜎𝑁,𝜎𝑆) is defined (see Figure 2. 18). The stress on a plane is presented in each point (𝜎𝑁,𝜎𝑆) in the stress space. Acceptable pairs (𝜎𝑁,𝜎𝑆 ) are given by the Equation (2.61) and all of these points have to fall inside a circle having center ( 1 2 (𝜎1 + 𝜎3), 0) and radius 1 2 (𝜎1 − 𝜎3), the large circle in Figure 2. 18. The two smaller circles are the points which lie outside the circles with center ( 1 2 (𝜎1 + 𝜎2), 0) and radius 1 2 (𝜎1 − 𝜎2). In conclusion, acceptable points in the stress space fall in the shaded region of Figure 2. 18. 41 Figure 2. 18: Mohr diagram in 2-D with appropriate points in stress space https://courses.eas.ualberta.ca/eas421/lecturepages/stress.html 2.3.3 P- and T-axis Using the method of first polarities of the P-arrivals for focal mechanism determination, four quadrants are created in which two are dilatational and two are compressional. The P- and T-axes lie in the center of the compressional and dilatational quadrants, respectively. The tension (T) and the pressure (P) axes are the axes of maximum extension (σ3) and the maximum compression (σ1), respectively. The i