Probabilistic seismic hazard analysis and site characterization of a power plant site in Chittagong, Bangladesh

A power plant is an essential structure that powers various sectors of the economy. It is vital to ensure the safety of such structures from natural hazards like earthquakes. The probabilistic seismic hazard analysis (PSHA) is the most advanced method of estimating earthquake hazards and risks. PSHA is the standard method for making site-specific, earthquake-resistant designs for important structures like a power plant. In the present study, the PSHA has been carried out for a power plant site in Chittagong using the latest available information on seismicity in the region. An updated earthquake catalog homogenized in uniform moment magnitude scale was prepared for the region. Seismicity analysis was carried out, and the seismicity parameters for the region were estimated from the frequency magnitude distribution plot. The peak ground acceleration (PGA) and the spectral acceleration at the bedrock for 10% and 2% probability of exceedance within 50 years were evaluated using different source models and attenuation relations in a logic tree framework. The site was characterized based on the soil investigation data from 31 boreholes, and amplification factors considering the local geology (SPT-N values) were obtained through ground response analysis. Subsequently, the surface level peak ground acceleration was estimated and depicted as hazard contour maps. The power plant site was observed to have a large spatial variation of seismic hazard at different locations.


Introduction
Chittagong is one of the major cities in Bangladesh, significantly impacting the country's economy.Since the southwestern region of Bangladesh lies in the active-tectonic region, it has experienced many large earthquakes in its history.A 225 MW combined cycle power plant was recently commissioned at Sikalbaha, Chittagong.A power plant is an essential structure that powers various sectors of the economy.Hence, its safety must be assured from natural hazards like earthquakes.To prevent structural damages caused by earthquakes, precautions must be taken to design structures to withstand earthquake forces estimated from seismic hazard studies.The probabilistic seismic hazard analysis (PSHA) is the most advanced method of determining the ground motion characteristics and their associated risk.According to US Nuclear Regulatory Commission (USNRC 2007), PSHA is the standard method for making site-specific, earthquake-resistant designs for the structures of a power plant.The uncertainties like the location, magnitude, and recurrence interval of an earthquake are considered in the methodology of PSHA.In the survey conducted by Goda et al. (2015) after the 2015 Nepal earthquake, it was observed that the local geology significantly impacted the building collapse pattern.In a particular area where alluvial and Pleistocene soil deposits are a few meters apart, it was found that the number of buildings that collapsed was twice in the alluvial deposit than in the Pleistocene deposit.Furthermore, liquefaction and differential settlement of structures were also detected close to this region.This indicates the amplification of the ground motion characteristics due to the local site conditions.Historically, similar observation has been made in the 1985 Mexico earthquake, the 1989 Loma Prieta earthquake, and the 2001 Bhuj earthquake, among many others.The liquefaction potential study of the Sikalbaha power plant site conducted by Sengupta and Kolathayar (2020) through a deterministic approach indicated a high susceptibility to liquefaction in the region.However, they used the amplification factor derived from broad NEHRP site classification.
The proximity to the convergent boundary of the Indian and Eurasian plates is one of the critical reasons for the high seismicity around the region (Morino et al. 2011(Morino et al. , 2014;;Steckler et al. 2018).There is also a significant prevalence of non-engineered construction, poor urban planning, and a high population density.These factors escalate the seismic hazard in the region drastically, creating a need for precise knowledge of the ground motion parameters.Several researchers have carried out seismic hazard studies for this purpose (Bhatia et al. 1999;Al-Hussaini and Al-noman 2010;Trianni et al. 2014, Al-hussaini et al. 2015;Das et al. 2016;Khan et al. 2018;Rahman et al. 2018).Bhatia et al. (1999) performed the probabilistic seismic hazard analysis of India and the adjoining regions (from 0 to 40°N and 65 to 100°E) under the Global Seismic Hazard Assessment Program (GSHAP).Al-Hussaini and Al-noman (2010) also performed the seismic hazard analysis with a probabilistic approach, employing several GMPEs.However, these GMPEs were not incorporated into a logic tree framework, leaving epistemic uncertainties possible.A probabilistic seismic hazard study to estimate the potential seismic risk along an oil pipeline in the Bay of Bengal was carried out by Trianni et al. (2014).Al-Hussaini et al. (2015) conducted a review study of Bangladesh's previous seismic hazard studies.They also performed the neo-deterministic seismic hazard analysis (NDSHA) and compared it with the previous SHA results.Das et al. (2016) estimated the probabilistic seismic hazard for a 100-year, 225-year, 475-year, 2475-year, and 10,000year return period for the northeast region of India, which also covered Bangladesh.Khan et al. (2018) performed the probabilistic seismic hazard assessment using a logic tree framework while incorporating the V s30 values over the entire country.Conversely, employing V s30 results for such a vast region produces conservative results due to the area's discontinuity and lack of sample data points.The seismic zonation map of Bangladesh, as illustrated by the Bangladesh National Building Code (BNBC 2017), is based on the maximum credible earthquake (MCE), with a 2475year return period, i.e., a 2% probability of exceedance in 50 years.Using this, design basis earthquake (DBE) can be found, which is 2/3 of the MCE and is equivalent to a ground motion with a 475-year return period or 10% probability of exceedance in 50 years.According to the standard, a PGA value of 0.19 g has been stipulated for Chittagong city.
There are several site response studies reported in the literature for different regions across the globe to estimate surface level seismic hazard, and the application of site response analysis is popular in recent studies due to its critical importance (Derras et al. 2020;Régnier et al. 2020;Adampira and Derakhshandi 2020;Pavel et al. 2020;Guzel et al. 2020;Kakhki et al. 2020).Koçkar and Akgün (2008) developed a geotechnical database for Ankara, Turkey, and integrated them with in situ characterization studies.Koçkar et al. (2010) and Koçkar and Akgün (2012) performed seismic site characterization of near-surface geologic materials for the Ankara basin and studied the site effects for the region.Yousefi-Bavil et al. (2022) proposed a 3D basin model to assess the site effects in the active tectonic regions of Gölyaka basin.
Raqib Al Mahmood and Hossain (2015) performed ground response analysis for Ganakbari, in Dhaka, Bangladesh using SPT-N values of nine boreholes.They used correlations to determine shear wave velocity values from SPT-N values.Kolathayar et al. (2013) evaluated the liquefaction potential of India based on SPT-N values required to prevent liquefaction using a probabilistic approach.They used NEHRP site classification based on satellite data to estimate site amplification.Vipin et al. (2013) conducted a performance-based study evaluating two important geotechnical effects of earthquakes, site amplification and liquefaction, for the state of Gujarat based on NEHRP site classes.Farazi et al. (2023) estimated shear wave velocity profiles for 19 seismic stations in Bangladesh using horizontal to vertical spectral ratio and reported that the region has deep soft soil deposits.
It can be observed that there is a lack of implementation of site-specific soil data in the evaluation of ground motion parameters in the Chittagong region of Bangladesh.Hence, an attempt has been made to evaluate the seismic hazard of a newly developed power plant site in Chittagong, employing a probabilistic approach and considering the local geology of the region.

Study area
Chittagong (Fig. 1) is Bangladesh's second most important city after Dhaka.It has military and economic importance due to the presence of the busiest ports in the country.In 2017, a 225 MW combined cycle power plant was commissioned in Sikalbaha, Chittagong.It lies along the banks of the Karnaphuli River in the north and is surrounded by other industrial establishments.The plant site lies within 22.323° N and 22.327° N latitudes and 91.863° E and 91.87° E longitudes.
According to the study by, Rahman and Siddiqua ( 2016), Chittagong primarily comprises sandstone, shale, and siltstones in its underlying geology.The city lies on the western edge of the Tripura-Chittagong folded belt, which is made up of several parallels and sub-parallel plunging folds (Fig. 2).
The major faults present in the city are the Sitakunda fault which cuts by the Karnaphuli River, the Tiger pass fault which runs along the NE-SW direction, and crosses the Karnaphuli fault almost orthogonally (Muminullah 1978).Khan and Chouhan (1996) suggested that the majority of the seismic events in the region are caused due to the subduction of the Bengal Basin region into the Burmese subplate.

Earthquake data processing
A detailed and comprehensive earthquake catalog is essential to any seismic hazard assessment.The earthquake catalog for this study was obtained from the US Geological Survey (USGS).Due to the non-uniformity of magnitude scales prevalent in such catalogs, the earthquake events were homogenized to the moment magnitude scale (M w ).
The homogenization process was performed using the relations proposed by Scordilis (2006), as mentioned below: The earthquake events of magnitude M w ≥ 4 in the past 120 years, i.e., 1900 to 2020, within a radius of 500 km from the site, are included in the catalog.Since the earthquake events are independent and obey the Poisson distribution, the insignificant foreshock and aftershock events must be removed.These extraneous events were removed by declustering using the Reasenberg (1985) algorithm incorporated in ZMap (Wiemer 2001).Eleven clusters (42 out of 1844 events) were found, and 1813 main events of M w ≥ 4 in the period were finally obtained from the raw data (Fig. 3).
The catalog's completeness was analyzed using the maximum likelihood method (Bender 1983;Aki 1965;Utsu 1999).This method uses earthquake events greater than the magnitude of completeness (Mc) to determine the "a" and "b" values (Gutenberg and Richter 1954).The magnitude of completeness is the minimum magnitude above which it can be reliably considered that the earthquake catalog is comprehensive.The b-value can be M w = 0.67 (±0.005)M s + 2.07 (±0.03) (for 3.0 ≤ Ms ≤ 6.1) M w = 0.99 (±0.02)M s + 0.08 (±0.13) (for 6.2 ≤ Ms ≤ 8.2) M w = 0.85 (±0.04) m b + 1.03 (±0.23) (for 3.5 ≤ mb ≤ 6.2) estimated by the relation given by the maximum likelihood method: where Δm bin is the bin size of the magnitude, M c is the magnitude of completeness, and the m mean is the average magnitude of the catalog.One of the most popular methods to estimate M c is the power-law curve fit of the frequency-magnitude distribution (FMD) plot, as Wiemer and Wyss (2000) suggested.The increment of the magnitude values or the bin size used in this study is 0.1.Figure 4 shows the frequency-magnitude distribution for the catalog along with the M c value.The seismicity parameters for the region were estimated as a = 8.657 and b = 1.21 ± 0.04.

Ground motion prediction equations (GMPE)
The ground motion models for this region are yet to be developed.Since the plant site lies in the active tectonic region, GMPEs developed for active tectonic areas in other countries have been used to evaluate seismic hazard.Attenuation models by Ghasemi et al. (2009)

Seismic hazard assessment
The probabilistic seismic hazard analysis was carried out as per the procedure developed by Cornell (1968).This study has incorporated two seismic source models: (a) linear source model and (b) gridded seismicity model.

Linear source model
The linear sources of India and the adjoining regions have been illustrated in the Seismotectonic Atlas published by the Geological Survey of India (GSI) (Dasgupta et al. 2000).SEISAT contains the details of faults and other linear sources in and around Bangladesh as well.It has been used as a credible source of information for faults, shear zones, lineaments, and geological features by various researchers in the past.Iyengar and Ghosh (2004) used it for their seismic hazard study of Delhi, Raghukant and Iyengar (2006) for Mumbai, Anbazhagan  2009) for South India.The map was first digitized, and the declustered events were layered on top of it in GIS software.The number of earthquake events and the maximum magnitude for each source were recorded.If an event is within 15 km of more than one source, then, the event is considered the source having the least distance to the epicenter.One significant limitation of the linear source model is that some of the events do not fall within the range of any fault.To overcome this shortcoming of the model, another seismic model proposed by Frankel (1995) is implemented in the logic tree.

Gridded seismicity model
The gridded seismicity model or the zoneless approach is used to overcome the limitations mentioned above of the linear source model.It also considers the Seismotectonic Atlas' inadequacy in the presence of hidden faults (Frankel 1995;Woo 1996;Martin et al. 2002).This method has been widely used in many research studies like Wahlström and Grünthal (2000), Lapajne et al. (2003), Jaiswal and Sinha (2007), Kalkan et al. (2009), Menon et al. (2010), Sitharam and Vipin (2011), Kolathayar andSitharam (2012), andSitharam et al. (2015).This model divides the entire area into several grids, and the number of events above a certain threshold magnitude is counted (M cut ).The recurrence interval of each magnitude range is noted, which is then smoothed using a Gaussian function to obtain the final activity of that grid cell.This study considered grids of 0.1° × 0.1° and threshold magnitude of M w = 4.0.The smoothening was done by the equation: where ni is the smoothed number of events in the i-th grid, n j is the number of events in the j-th grid, ∆ ij is the distance between the i-th and j-th cell, and c is the correlation distance to be considered for location uncertainties.

Logic tree structure
To address the issue of epistemic uncertainty, the logic tree framework was implemented to minimize errors (Bommer et al. 2005;Budnitz et al. 1997;Stepp et al. 2001).A logic tree is a branched structure that presents different models at each node.Appropriate weights are allocated at each node, with the condition that the sum of all the weights at a node should be unity.In the present study, the main factors of the logic tree are (a) different seismic sources, (b) evaluation of the maximum earthquake magnitude, and (c) chosen GMPEs in the study.
In this investigation, both the seismicity models are equally relevant; hence, a weightage of 0.5 has been assigned to both.The maximum earthquake magnitude is calculated using two approaches: for the linear source model, (Akkar and Bommer 2010) the actual observed maximum magnitude and (Al-Hussaini and Al-Noman 2010) the actual maximum magnitude plus 0.5, whereas for the gridded seismicity model, (Akkar and Bommer 2010) the actual observed maximum magnitude and (Al-Hussaini and Al-Noman 2010) the maximum magnitude obtained from the process described by Kijko (2004).Both these approaches were given equal weights in the logic tree.Figure 5 shows the final logic tree structure.

Local site effects
Local geology significantly impacts the ground motion characteristics obtained at the bedrock level.Since the area of interest in the study is quite small, it was possible to  for evaluating the surface level PGA.The SPT values of all the boreholes were given as input in the DEEPSOIL software v7 (Hashash et al. 2016), from which the amplification factors for each borehole were obtained.The surface PGA was evaluated by multiplying the bedrock PGA with the amplification factors.Figure 6 shows the plant area and the points where standard penetration measurements were taken.Figure 7 summarizes SPT-N values for BH 119 to BH 133. Figure 8 shows a typical bore log for borehole 125 and also the corresponding variation of SPT-N value with depth.The site's geotechnical data, including SPT-N values, were available for different boreholes at various depths.However, shear wave velocity values were not available from the field tests.Hence, in the present study, the established correlations were to estimate shear wave velocity values at different depths for different boreholes.The correlations were chosen based on the guidelines for estimation of shear wave velocity profiles, published by Pacific Earthquake Engineering Research Center (Wair et al. 2012).Two well-established correlations were employed to estimate V s30 from SPT-N values.The V s30 values for each borehole were calculated from the available SPT-N values by the conversion relations given by Seed et al. (1986) and Yoshida et al. (1988).The results from both regions indicate that the variation is marginal.A slightly lower value of shear wave velocity was obtained from the relation proposed by Seed et al. (1986), and hence, the same was considered for site response analysis, to be on the conservative side.An artificial acceleration time history was generated for the power plant site considering the estimated PGA and target response spectrum using the software SeismoArtif which can simulate ground motions compatible with given response spectra (Gasparini 1976).The site-specific input motion generated is shown in Fig. 9.In the DEEPSOIL program, this site-specific synthetic accelerogram developed for the power plant site was given as the input motion.The software then plots several ground motion characteristics and in situ soil properties.It also generates the response spectrum at the bedrock and the surface level.The PGA value of the surface is then divided by the PGA at the bedrock to get the amplification factor for that specific borehole location.The ground response analyses were carried out for all the boreholes considering their soil properties.Sample analysis for a typical borehole (BH-125) is presented here.Figure 10 shows the ground motion characteristics and in situ soil properties generated by DEEPSOIL at BH-125. Figure 11 illustrates the response spectrum at bedrock (input motion) and surface level (Layer 1) due to the given input motion at BH-125.The amplification factors for all borehole locations were calculated through a similar analysis.
The above procedure was repeated for all the boreholes to develop a site amplification contour.Table 1 illustrates the location of each borehole in the site, along with their amplification factors.

Results and discussion
The probabilistic seismic hazard analysis for the power plant site was carried out using appropriate source models and GMPEs in a logic tree framework.Hazard contour maps Page 13 of 15 542 about 0.19 g for 0.1 s period and 0.04 g for 1 s period can be expected.In the case of 2475-year return period, a value of 0.29 g for 0.1 s period and 0.06 g for 1 s period can be anticipated.From the hazard maps, a general observation can be made that the bedrock PGA is relatively uniform over the entire power plant site area.Table 2 compares the results of this study with previous investigations and available codal provisions.From the table, it can be established that the results of this study are reliable.
Figure 18 shows the variation of the amplification factor over the entire site, which was obtained from the results of the DEEPSOIL analyses of 31 borehole sites (Table 1).The surface level PGA values at different locations were obtained for 475 and 2475 years return periods (Figs.19 and 20) by multiplying the bedrock motion with the amplification factor at corresponding locations.On comparison of the obtained PGA values with those stipulated by the BNBC-2017 (0.19 g), it can be concluded that the PGA values specified by the code match reasonably well.However, a site-specific study is always advisable for essential structures, such as a power plant, since the local geology significantly impacts the strong ground motion.The surface level hazard was observed to vary across the power plant site (Figs. 19 and 20), whereas the bedrock level hazard was relatively uniform.This indicates the importance of local site effects on seismic hazards.

Conclusion
The probabilistic seismic hazard analysis of a power plant site at Sikalbaha, Chittagong, was carried out using the Cornell-McGuire approach using multiple seismic source models and GMPEs.The study was carried out with the most recent data, and a tool such as a logic tree was implemented to rule out epistemic uncertainties.It was found that a bedrock PGA of 0.11 g can be expected for a return period of 475 years and 0.15 g for a return period of 2475 years at the power plant site.The results were in agreement with previous studies and recent codal provisions.
Based on the soil investigation data available from 31 boreholes, detailed ground response analyses were carried out for all the borehole locations.It was found that the spatial variation of the soil profile was significant, and hence, the amplification factors too varied across the site.The estimated hazard at the bedrock level was brought to the ground level using the amplification factors for different locations at the site.The surface level PGA for return periods of 475 years and 2475 years were as high as 0.12 g and 0.18 g, respectively, in the northeastern part of the site.Unlike the variation of bedrock level PGA, which had a uniform spatial variation, the surface level PGA varies significantly at different locations.This indicates that the effect of amplification is profound at the site.Extensive geological information on the seismic sources will aid in better estimation of hazards.Development of regionspecific ground motion prediction equations using actual strong motion data for the region will go a long way in more reliable seismic hazard estimation.Geophysical surveys to estimate shear wave velocities can give a more reliable estimate of the site amplification than geotechnical data correlations.

Fig. 1
Fig. 1 Map showing the city of Chittagong with important locations

Fig. 2
Fig.2Major faults and geological features in and around Bangladesh(Khan and Chouhan 1996)

Fig. 4
Fig. 4 Frequency-magnitude distribution plot depicting the magnitude of completeness Fig. 5 Logic tree framework

Fig. 6 Fig. 7 Fig. 8
Fig. 6 Power plant area (unshaded area) along with points where SPT tests were conducted

Fig. 9
Fig.9The artificially generated acceleration/velocity/displacement time history for the power plant site

Fig. 10
Fig. 10 Ground motion characteristics and in situ soil properties generated by DEEPSOIL at BH-125

Fig. 11
Fig. 11 Response spectrum at bedrock (input motion) and surface level (Layer 1) due to the site-specific input motion at BH-125

Fig. 15
Fig. 15 Sa contour at 1 s for a 475-year return period

Table 1
Location of each borehole on the site, along with the evaluated amplification factors