Probabilistic seismic hazard assessments for Myanmar and its metropolitan areas

Although Myanmar is an earthquake-prone country, there has not been proposed an official national seismic hazard map. Thus, this study conducted a probabilistic seismic hazard assessment for Myanmar and some of its metropolitan areas. Performing this assessment required a set of databases that incorporates both earthquake catalogs and fault parameters. We obtained seismic parameters from the International Seismological Centre, and the fault database includes fault parameters from paper reviews and the database. Based on seismic activities, we considered three categories of seismogenic sources—active fault source, shallow area source, and subduction zone source. We evaluated seismic activity of each source based on the earthquake catalogs and fault parameters. Evaluating the ground-shaking behaviors for Myanmar requires evaluation of ground-shaking attenuation; therefore, we validated existing ground motion prediction equations (GMPEs) by comparing instrumental observations and felt intensities for recent earthquakes. We then incorporated the best fitting GMPEs into our seismic hazard assessments. By incorporating the Vs30 (the average shear velocity down to 30 m depth) map from an analysis of topographic slope, we utilized site effect and assessed national probabilistic seismic hazards for Myanmar. The assessment shows highest seismic hazard levels near those faults with high slip rates, including the Sagaing Fault and along the Western Coast of Myanmar. We also assessed seismic hazard for some metropolitan cities, including Bagan, Bago, Mandalay, Sagaing, Taungoo and Yangon, in the forms of hazard curves and disaggregation by implementing detailed Vs30 maps from micro-tremor surveys. The city-scale assessments show higher hazards for sites close to an active fault or/and with a low Vs30, demonstrating the importance of investigating site conditions. The outcomes of this study will be beneficial to urban planning on a city scale and building code legislation on a national scale.


Introduction
Due to the collision between the Eurasian and Indo-Australian plates, Myanmar is an area with a high seismic hazard level.According to the earthquake parameters summarized by the International Seismological Centre (ISC, http:// www.isc.ac.uk, shown in Fig. 1), around 140 events with M ≥ 3.0 have taken place in Myanmar and its vicinity every year from 1990 to 2019.In addition, some significant events with M ≥ 6.0 resulted in severe casualties and building damage, such as the 2011 M w 6.9 earthquake in the State of Shan and the 2012 M w 6.1 earthquake in Thabeikkyin.As a result, it is desired to understand the seismic hazard potential for Myanmar, especially for some populous cities, for which a probabilistic seismic hazard assessment would be a practical approach.
A probabilistic seismic hazard assessment (PSHA; e.g., Cornell 1968) illustrates the probability of exceedance of ground motion in a specified period.In other words, seismic hazard in the form of strong ground shaking is quantified by summarizing the impacts from all the seismogenic sources near the site(s) of interest.PSHAs have been widely applied to engineering purposes, for example, infrastructure sites selection, safety evaluation for nuclear power plants, building code legislation, and insurance premium determination.Thus, many PSHAs have been proposed for various spatial scales.For example, the Global Earthquake Model (https:// www.globa lquak emodel.org/ gem) has proposed a global seismic hazard map by homogenizing seismic hazard levels across the seismic models from various studies (Pagani et al. 2020); Ornthammarath et al. (2020) and Peterson et al. (2007) have proposed a regional seismic hazard assessment for southeast Asia.Somsa-Ard and Pailoplee (2013) proposed both deterministic and probabilistic seismic hazard assessments considering seismic activity determined by earthquake catalog.Thant (2014) conducted a PSHA for the Yangon region, which can be regarded as one of the first assessments in a cityscale for Myanmar metropolitan.Although the assessments conducted in these studies include Myanmar, these assessments lacked detailed information, such as active fault parameters, strong motion attenuation behaviors and site amplification factors, on this country.Pailoplee et al. (2009) proposed a PSHA based on earthquake catalog.Afterward, earthquake catalog has been updated, some significant events took place in Myanmar (Fig. 1), and seismicity activity has been investigated by previous studies (e.g., Wang et al. 2014), providing good constraint on subsequent hazard assessments.
Conducting a reliable PSHA requires understanding seismogenic sources, ground motion attenuation behavior, and site amplification of target sites.To propose a reliable seismic model for each seismogenic source, we evaluated the possible recurrence interval of each source by fitting with the observations from a database that comprises active fault and earthquake parameters.To accommodate the path effect for various tectonic regimes, we incorporated proper ground motion prediction equations (GMPEs) selected by a series of tests comparing strong ground motion observations of some representative earthquakes.To incorporate site amplification into the PSHA, we introduced the average shear-wave velocity for the upper-most 30 m depth (V s 30 ) obtained by topographic slope on a regional scale and micro-tremor survey on a city scale.Considering the above-mentioned procedures, we assessed seismic hazard for Myanmar and some of its metropolitan areas in the forms of a hazard map, hazard curves, and disaggregation.We discussed the outcomes of the PSHA and potential applications.

Probabilistic seismic hazard assessment
A traditional PSHA was first developed by Cornell (1968), assuming that earthquakes are independent of each other and earthquake probability P follows a Poisson process, expressed as where t is the period of interest; and ν is the average return period in a seismogenic source, usually in the form of (1) where N s i=1 is the summation of the probability of N earthquake events; P[Y > y * [m, r]] is the probability of ground-shaking; Y exceeds y * for earthquakes of magni- tude m at distance r; f M i and f R i are probability density functions of magnitude and distance, respectively, and N i is the annual seismic rate.For example, the hazard level for probabilities, p(t, ν) , of 2% and 10% in a period, t, of 50 years correspond to return periods, ν , of 2,475 and 475 years, respectively.The hazard levels in these two return periods are most common standards for the design codes of infrastructures and common buildings, respectively.To describe N i , we applied the trun- cation exponential model and the characteristic model, described in the following sections.For PSHA calculation, we implemented the OpenQuake-engine hazard calculation software (Pagani et al. 2014).

The truncated exponential model
The truncated exponential model is based on the Gutenberg-Richter law (Gutenberg and Richter 1944), which describes magnitude-frequency distribution as an inverse power law.By further considering the upper and lower bounds of magnitudes, the truncated exponential model can be expressed as where Ṅ is the seismic rate with magnitude larger or equal to M , and a and b are constants obtained by regres- sion according to seismic activity or estimation according to fault parameters (Fig. 2).

The characteristic model
While the truncated exponential model suggests seismicity activity in a seismogenic source within a magnitude range, previous studies (e.g., Wang et al. 2014;Chan et al. 2020) inferred the seismic hazards induced by specific faults using the seismicity rates of the characteristic earthquake with the maximum magnitudes.
The occurrence rate of the characteristic earthquake, Ṅ , could be estimated with the interseismic slip rate of a fault ( Ṡ ) and displacement of a characteristic earthquake (D):  S1) based on the exponential model

Seismic source with modelled seismic activity
In this study, the seismogenic sources were divided into three categories: active fault sources, shallow area sources, and subduction interface sources.These are detailed in the following sections.

Active fault sources
To propose a reliable seismic model for active fault sources, we accessed the fault database summarized by Earth Observatory of Singapore (EOS).This fault database was originally proposed by Wang et al. (2014) and subsequently incorporated active faults mapped by remote sensing in northeastern India, Bangladesh, Thailand, Laos, Vietnam and Southern China (Fig. 3).The fault parameters of each active fault were shown in Additional file 1: Table S1.Based on the active fault map, the studied region could be separated into different tectonic regimes.The Sagaing Fault with a high slip rate runs across the center of Myanmar and can be associated with several significant events that took place in Myanmar (Fig. 1).There are many faults with northeast-southwest and northwest-southeast trending with the potential of earthquakes with M > 6.0 in the Shan Pleateau adjacent to China, which can be associated with the clockwise rotation of the Himalayan tectonic system (Shi et al. 2018).
The Three Pagodas fault system (Fault 252 and 254) along the border of Myanmar and Thailand could produce earthquakes with a magnitude higher than 7.0.
Based on the active database, the parameters of 430 active faults and the surface alignments were obtained (Fig. 3).To model the recurrence interval of each active fault, we implemented both the truncated exponential model and the characteristic model.
For our truncated exponential model, we assumed a minimum magnitude of M6.5 for all active faults.Since the majority of faults in Myanmar are strike-slip faults, we believe the surface rupture of the fault can be identified.In addition, we determined the maximum magnitude according to fault length and the empirical formula proposed by Wells and Coppersmith (1994).
We used the truncated exponential model to derive the occurrence rate for each magnitude bin: The ratio for the recurrence rate for magnitude difference of 0.1 can be expressed as Here, we assumed a fixed b value of 1.0 for this model and used the scaling law proposed by Wells and Coppersmith (1994) to obtain displacement of an earthquake D(M) with a magnitude M: The c and d for different focal types are shown in Table 1.Considering Eqs. 6 and 7, we obtained: (5) where Ṡ is the sum of the slip rate contributions from all the magnitude bins (shown as the color for each fault in Fig. 3).

Fig. 4 continued
Take Fault 117 (a segment of the Sagaing Fault, shown in Fig. 3) as an example (Fig. 2), the interseismic slip rate of this fault segment is 18 mm per year.We assumed the magnitude range is between 6.5 and 8.1.Since the Sagaing Fault is of a strike-slip type, we implemented − 6.32 and 0.90 as the values of c and d, respectively (Table 1).Considering all the parameters, we then obtained the seismic rate for each magnitude bin for the fault.
In addition to the exponential model, the characteristic model was implemented for active fault sources.The characteristic displacement of an earthquake with a magnitude M for each fault was obtained using the scaling law (Eq.7).By considering the slip rate of each fault from our database, the occurrence rate was obtained using Eq. 4.
Our PSHA considered both the exponential and the characteristic models and assumed equal weighting (50% for each) in each logic tree branch for the assessment.

Shallow area sources
In addition to the active fault sources, some of the crustal earthquakes could not be associated with any identified fault.We thus introduced shallow area sources to model the seismic activities.We divided the study region into 14 area sources (Fig. 3) by tectonic setting and seismicity activity (Wang et al. 2014).Zone 1 represents the Sagaing fault system, a dextral slip associated with the northward translation of India along a 1400 km span centered on Myanmar (e.g., Curray et al., 2002;Vigny et al. 2003 (Wang et al. 2014); and Zones 10-13 represent the Shan domain with regional bookshelf faulting (Shi et al. 2018) between the Sagaing (Faults 112-126) and the Red River fault (Fault 10).
To model the seismic rate for each area, we accessed the ISC earthquake catalogue (http:// www.isc.ac.uk) to obtain earthquake parameters, including occurrence time, longitude, latitude, depth, and magnitude, for the earthquakes since 1904.According to the earthquake epicenters, the corresponding seismicity in each zone was sorted.To follow a Poissonian process (Eq.1), earthquakes in the catalog were assumed to be independent of each other.Thus, we implemented the declustering approach proposed by Gardner and Knopoff (1974) to remove foreshocks and aftershocks from a catalog.Based on the declustered catalog, we fit the observations through the truncated exponential model (Eq.3) to describe the seismicity rates of the shallow area sources (Table 2).
Take Zone 1 as an example, we first indicated complete part of the catalog for both in time and magnitude.In the earthquake catalog since 2003, we obtained magnitude of completeness (M c ) of 3.5 based on the maximum curvature method (Wiemer and Wyss 2000) according to distribution of the annual seismicity rate as a function of magnitude (Fig. 4a).Considering the complete part of the catalog, we fit the seismicity rate for magnitudes between 3.5 and 5.5 through the regression of a truncated model line and obtained a and b values of 3.63 and 0.87, respectively (Fig. 4a).The completeness and maximum magnitude implemented for regression of a truncated model line were summarized in Table 2. Based on this procedure, we implemented to all of the zone to determined M c , catalog completion year, a and b values (Fig. 4).Considering minimal magnitude for active fault sources is 6.5, we implemented maximal magnitude of 6.4 for shallow area sources to illustrate seismic activity for all of magnitude range.

Subduction sources
To the western offshore area of Myanmar, the Indian Plate is subducting northeastward beneath the Eurasian Plate.The Arakan thrust and the Andaman trench are convergent plate boundaries with the potential of megathrust earthquakes, for example, the 2004 M9.1 Sumatra-Andaman earthquake (Ammon et al. 2005).To model the seismicity along the subduction interface, referred to as the Arakan megathrust (illustrated in Fig. 3), we considered Faults 144, 175, and 176 with short recurrence intervals for earthquakes with magnitude ≥ 7.5.To model the seismicity rates for the three subduction interface sources, we followed the procedure for the active fault source to consider both the exponential and the characteristic models and assumed equal weighting (50% for each) in a logic tree.
Besides earthquakes along the megathrust, we considered seismicity activity from three area sources (Zones 3-5, shown in Fig. 3), with Zones 3 and 4 representing the area between the Arakan trench and the north of the Arakan range and Zone 5 representing the area along the Andaman Trench.We proposed the truncated exponential model (Eq.3) to describe the seismicity rates of the three area sources (detailed in Table 2).

Ground motion prediction equations
In addition to the seismic source model, another crucial factor for seismic hazard assessment is the proper evaluation of ground-shaking attenuation.To do this in different tectonic regimes, we conducted a series of tests that constrained estimated ground motions by actual felt intensities of earthquakes.To model scenario groundshaking levels, we used GMPEs, which estimate the ground shaking attenuation as a function of source characteristics and path and site effects.We tested the credibility of each GMPE from the Ground-Shaking Intensity Models (GSIM, https:// docs.openq uake.org/ oq-engine/ master/ openq uake.hazar dlib.gsim.html) database maintained by the Global Earthquake Model (GEM, http:// www.globa lquak emodel.org/).To comprehend different attenuation behaviors between subduction and crustal events, we conducted three sets of GMPE tests.For crustal earthquake sources, we implemented intensity reports summarized by "Did You Feel It?" (http:// earth quake.usgs.gov/ earth quakes/ dyfi/) for the 2011 M6.8 Tarlay and 2012 M6.8 Thabeikkyin earthquakes (shown in Fig. 5a, b, respectively).Considering earthquake parameters (i.e., earthquake magnitude, rupture length, shown in Fig. 5) and various GMPEs, we estimated possible ground shaking intensity (converted from PGA) at each site with intensity reports and reported the averaged differences with observations.An ideal case presents that estimated intensity perfectly match observations.In the both crustal cases, we concluded the ground-shaking behaviors of the crustal earthquake best fit with the GMPEs proposed by Akkar and Cagnan (2010), shown as where c 1 is a constant as a reference magnitude, fixed to be 6.5, F N and F R are variables for the influence of faulting, taking values of 1 for normal and reverse faults and zero otherwise.The coefficients of the regressions are given in Additional file 1: Table S2. (9a) For subduction events, we compared macroseismic observations for the 1975 M6.5 Bagan earthquake with a depth of 84 km summarized by Lin Thu Aung (personal communication).We concluded the groundshaking observations (Fig. 5c) best fit with the GMPEs proposed Atkinson and Boore (2003), shown as where Y is peak ground acceleration or 5% damped pseudo-acceleration in cm/sec random horizontal component, M is moment magnitude, fn( ) with D fault being the closest distance to fault surface, in kilo- meters (use h = 100 km for events with depth > 100 km) and is a near-source saturation term, assumed to be = 0.00724 × 10 0.507M and ( 10) g = 10 (1.2−0.18M) for interface events, g = 10 (0.301−0.01M) for in-slab events.sl = 1.for PGA rx ≤ 100cm/sec 2 or frequencies ≤ 1Hz sl = 1.− (f − 1)(PGA rx − 100.)/400.for 100 sl = 1.− (PGA rx − 100.)/400.for 100 < PGA rx < 500cm/sec 2 (f ≥ 2HzandPGA) sl = 0. for PGA rx ≥ 500cm/sec 2 (f ≥ 2HzandPGA), PGA rx is predicted PGA on rock (NEHRPB) , in cm/sec and σ is standard deviation of residuals, equal to σ 2 1 + σ 2 2 where 1, 2 denote estimated intra-and interevent variability, respectively.The coefficients of the regression are given in Additional file 1: Table S3.
We implemented the best fitting GMPEs to our PSHA.

Site effect
A site with soft soil condition could result in amplification of strong ground shaking (Lermo and Chávez-García, 1993).For engineering applications, a seismic hazard assessment usually implements a GMPE with a term of V s 30 , that is, a lower V s 30 (i.e., a site with soft soil at subsurface) infers a larger site amplification, and vice versa.We introduced two sets of V s 30 databases at different spatial scales.We derived a regional V s 30 map (shown in Fig. 6) from empirical relationships from the Shuttle Radar Topography Mission digital elevation models (Allen and Wald 2009).In addition, to better assess seismic hazard on the city scale, the V s 30 data for the five metropolitan areas were obtained by the micro-tremor array measurements.These detailed distributions of V s 30 provide a better constraint on the site conditions (Fig. 7).In total, we obtained 408, 660, 504, 302, and 83 sites for micro-tremor measurements for Sagaing, Mandalay, Taungoo, Bago, and Yangon metropolitan areas, respectively (shown as grey dots in Fig. 7).
Taking the V s 30 map in Sagaing (Fig. 7a) as an example, the V s 30 on the Sagaing Hill (V s 30 = 1230 m/sec) is significantly higher than other sites, inferring an insignificant site effect.Downtown Sagaing is in the region with a low V s 30 of 320 m/sec, inferring the potential amplification in seismic hazard of the site.

Probability seismic hazard assessment for Myanmar and its five metropolitan cities
Summarizing the source model, GMPE, and site effect discussed in the preceding section, we assessed the seismic hazard for Myanmar and some of its metropolitan areas in the form of a seismic hazard map, hazard curves, and disaggregation.

Hazard maps
To demonstrate the spatial distribution of seismic hazard in Myanmar, we presented seismic hazard maps for Myanmar and its vicinity (Fig. 8) in the forms of 2% and 10% probabilities of exceedance in 50 years, corresponding to the recurrence intervals of 2475 and 475 years, respectively (according to Eq. 1).The hazard maps suggest that the regions close to the faults with high slip rates or short recurrence intervals have high seismic hazard levels, for example, along the Sagaing Fault across The Sagaing Fault passes through Sagaing and Bago and near to Mandalay and Taungoo, which brings the level of peak ground acceleration (PGA) higher than 1.0 g for a recurrence interval of 475 years.The arc-shaped subduction system extending from the southwest coast of Myanmar to the northwest also has a high seismic hazard, though its PGA intensity is lower than that of the Sagaing Fault for the same recurrence interval.Since the surface alignment of the subduction zone (i.e., trench) took place in the offshore region, it is necessary to be aware of the potential of tsunamis.
To evaluate site amplification, we compared the hazard level differences considering the V s 30 map of Wald and Allen (2009) and a fixed V s 30 of 750 m/sec as a site on engineering bedrock (Fig. 9).Comparing the distribution of V s 30 (Fig. 6), the lower V s 30 is, the higher the seismic hazard will be.There are several metropolitan areas located in regions with a significant site effect, for example, Mandalay, Sagaing, Taungoo, Bago, and Yangon; thus, it is important to assess in detail the seismic hazards on a city scale.

Hazard curves
Since several populous metropolitan areas are in regions of high seismic hazard due to site effect (open circles in Fig. 9), we conducted detailed seismic assessments of them.For these analyses, we included V s 30 maps obtained from micro-tremor-array measurements (Fig. 7).Our seismic hazard estimations based on these data appear in the form of hazard curves (probability as a function of ground-shaking level) for representative sites in each city (Fig. 10).Higher hazard levels for some sites correlate with low V s 30 , for example, Dagon University in Yangon (Fig. 7e) and Downtown Sagaing City (Fig. 7a) and/or their proximity to active seismogenic structures, for example, Shwemawdaw Pagoda next to the Sagaing Fault in Bago (Fig. 7d).
According to the hazard curves in the form of PGAs for five metropolitan areas (Fig. 10a), high hazard levels were expected in Sagaing, Mandalay, and Bago with a 2% probability of exceedance in 50 years higher than 1.0 g.The high hazard level in the three cities could be associated with their proximity to the Sagaing Fault with a large magnitude and high seismicity rate (Fig. 2).By contrast, Yangon has lower hazards among the five metropolitan areas due to its remoteness from active faults (Fig. 3).
To assess seismic hazard for the potential risk to buildings, in addition to PGA, we demonstrated seismic hazard in pseudo-spectral acceleration (SA) for various periods (damping ratio = 5%).Spectra acceleration at 0.3 s, denoted as "SA(0.3)",has a natural period similar to that of low-storied buildings while high-storied buildings for SA at 1.0 s [SA(1.0)].Through SA analysis, we can evaluate the seismic risks of buildings with different heights.According to the hazard curves for SA(0.3)(shown in Fig. 11b), significantly high hazard levels exist for most of the metropolitan areas (except

Fig. 1
Fig. 1 Distribution of seismicity in Myanmar and its vicinity.Stars and dots represent earthquakes with Mw ≥ 6.0 and Mw ≥ 3.0, respectively, summarized in the ISC catalog from 1990 to 2019

Fig. 2
Fig. 2 Seismicity rate model for different magnitudes for Fault 117 (central segment of the Sagaing Fault, corresponding alignment shown in Fig. 3 and parameters shown in Additional file 1: TableS1) based on the exponential model

Fig. 3
Fig. 3 Distribution of area and active fault sources in Myanmar and its vicinity.The color of each fault denotes its slip rate.The modeled seismicity rate for each area source is represented as a and b values of the Gutenberg-Richter law in Table2.Fault parameters of each active fault are presented in Additional file 1: TableS1

Fig. 4
Fig.4Seismicity rate model for each zone (the polygon of each zone is shown in Fig.3).Based on the complete part of the catalog in each region, the model is obtained by the regression of the Gutenberg-Richter law (red solid line) with standard deviations (red-dashed line) ); Zones 3-5 represent the Sunda megathrust at the Dhaka, Ramree, and Coco-Delta segments, respectively, defined by Wang et al. (2014); Zone 8 represents seismicity between the Himalayan Frontal Thrust (Fault 155) and the Dauki Fault (Fault 145, the boundary of the Shillong Plateau, Morino et al. 2011); Zone 9 represents the Naga domain, east of the Shillong plateau and the northern edge of the Burma plate

Fig. 5
Fig. 5 Observed and predicted strong ground motions for (a) the 2011 Tarlay, (b) 2012 Thabeikkyin, and (c) the 1975 Bagan earthquakes.The earthquake parameters and implemented GMPEs were denoted in each case

Fig. 6
Fig. 6 Distribution V s 30 in the study region, estimated by topographic slope (Wald and Allen 2009).The locations of some metropolitan areas were indicated

Fig. 8
Fig. 8 Hazard maps of Myanmar and its vicinity for (a) 2% and (b) 10% probabilities of exceedance in 50 years (in PGA, in g).Site amplification was considered for the assessment based on the V s 30 map shown in Fig. 6

Fig. 9
Fig.9Difference in hazard levels considering the V s 30 map ofWald and Allen (2009) and a fixed V s 30 of 750 m/sec as a site on engineering bedrock.The hazard is presented in 2% probability of exceedance in 50 years (in PGA, in g)

Fig. 11
Fig. 11 Hazard curves in (a) PGA, (b) SA(0.3), and (c) SA(1.0) in 50 years for the five metropolitan areas.We presented the hazard at the site with the highest hazard level at each metropolitan area-Dagon University for Yangon, Amarapura for Mandalay, Downtown Bago for Bago, Taungoo government office for Taungoo, and Downtown Sagaing for Sagaing

Fig. 12
Fig. 12 Disaggregation analysis of PGA with 10% probability of exceedance in 50 years for (a) Sagaing, (b) Mandalay, (c) Taungoo, (d) Bago, and (e) Yangon.We presented the hazard at the site with the highest hazard level at each metropolitan area

Table 1
Wells and Coppersmith (1994)lations between slip rate and magnitude proposed byWells and Coppersmith (1994)

Table 2
Parameters of each area source, including source category, a and b values, magnitude of completeness (M c ), maximum magnitude (M max ) and catalog completion year Note that M c and M max are implemented for regression of a truncated model line.The maximal magnitude of all zones for PSHA is fixed to be 6.5 ASC Active Shallow Crust source, SZ Subduction Zone source