Probabilistic seismic hazard assessment for offshore structures in Andaman Sea

: A set of probabilistic seismic hazard maps for Offshore structures in Andaman sea has been derived using procedures developed for the latest US National Seismic Hazard Maps. In contrast to earlier hazard maps for this region, which are mostly computed using delineated seismic source zone, the presented maps are based on the combination of smoothed gridded seismicity, crustal-fault, and subduction source models. The ground motion hazard map is presented over a 10 km grid in terms of peak ground acceleration and spectral acceleration at 1.0 undamped natural periods and a 5% critical damping ratio for 10 and 2% probabilities of exceedance in 50 years, which have generally been used for Seismic Analysis and Design of Offshore Structures.


INTRODUCTION
The Andaman Sea is situated in an active backarc basin lying above and behind the Sunda subduction zone, which is the Indo-Australian and Eurasian boundary zone comprise the convergent margins, including the Burma oblique subduction zone, Andaman thrust and Sunda arc, to the North West, west and south, respectively. Within Andaman Sea, several earthquakes with magnitude greater than 4 have been observed during the years 1964 to 2017. However, the earthquake activity rate is much lower than those occurred near plate boundary. For this study, the southern part of Andaman Sea (Blue dash line and the study area is bounded latitude 5° to 10° and longitude 94° to 98°, Figure 1) is of special interested due to ongoing human activities for gas pipeline and offshore facilities, several platforms and subsea gas pipelines have been developed offshore. In addition, onshore supports such as control and maintenance centre along with gas metering stations have been developed in this region. The largest instrumental earthquake within this zone is Mw 6.3 on 16 May 1933 at 10 km depth. In addition, this region is situated about 200 and 300 km from the Sumatran faults and Sumatra subduction zone, respectively. And these seismic active structures have historically produced moderate to large earthquakes with long-period ground motions that were felt in highrise buildings in Singapore and Malaysia (Pan and Sun, 1996;Pan, 1997;Pan et al., 2001). A probabilistic seismic hazard assessment for offshore structures located in this area has been carried out. ground motion parameters include: horizontal Peak Ground Acceleration (PGA) and Spectral Acceleration (SA) values at the project site with the return periods of 475 and 2475 years.

SEISMOTECTONIC SETTINGS
The Indo-Australian and Eurasian boundary zone comprise the convergent margins, including the Burma oblique subduction zone, Andaman thrust and Sunda arc, to the North West, west and south, respectively. The plate kinematics of the Sumatran region is, in a broad sense, the simple interaction of the Indian-Australian and Eurasian plates (Figure 1). However, in detail it is much more complex than that. Deformation rates across these plate boundaries are variable. The observed seismicity and seismotectonic settings of these plate boundaries clearly indicate the capability of producing large events, where the 26 December 2004 earthquake occurred. A convergence rate of 65-70 mm/year as a result of Australia moving toward South East Asia is reported by McCaffrey (1996).
Deformation of the overriding plate leads to larger complexities in plate motions. Sumatra sits at the southwestern edge of the Sunda plate (Bird 2003), which moves at a few millimeters per year to a centimeter per year eastward relative to Eurasia (Chamot-Rooke & Le Pichon 1999, Bock et al. 2003 (Figure 1). The resulting convergence between the Sunda plate and the oceanic plates to the southwest is somewhat slower than it would be relative to Eurasia. The rate and direction of subduction of the lithosphere under the Sunda forearc, however, are further modified by the independent motion of the forearc. Fitch (1972) explained the presence of the Sumatran fault and other similar faults inboard subduction zones by the process now known as slip partitioning. That is, in some cases of oblique subduction where the two plates do not converge at a right angle to the strike of the trench, it requires smaller overall shear force to share the shearing (trench-parallel) component of the relative motion between two separate faults instead of on one fault. In the case of partitioning, one fault is the subduction thrust, which takes up all of the trench-normal slip (the dip-slip component) and some fraction of the trench-parallel slip (the strike-slip component). A second fault, within the overriding plate and commonly strike-slip in nature, takes up a portion of the trench-parallel motion. The subduction thrust and strike-slip fault isolate a wedge of forearc called the sliver plate. The slip rates on the separate faults can be inferred from their geometries and knowledge of the overall convergence.
The motion of the Sunda forearc (sliver plate) is not known well, particularly in the Andaman section, and hence the subduction vector is highly uncertain. Earlier estimates of the relative motions assuming a rigid forearc sliver plate failed to predict convergence in the Andaman section of the trench, which probably indicates, as is now accepted, that there is extensive internal deformation within the forearc. One possibility, evident in the change in earthquake slip directions along the margin, is that the Andaman section (called the Burma plate) and Sumatran sections of the forearc move independently with a break near 6° to 7°N (Subarya et al. 2006).

EARTHQUAKE CATALOGUE
The earthquake catalogue for current study is composed of instrumental earthquake records from different international earthquake observatories including:

International Seismological Centre-Global
Earthquake Model (ISC-GEM) (1907-2009), 2. Engdahl (EHB)'s earthquake catalog , 3. USGS/NEIC preliminary earthquake database , and 4. Global centroid moment tensor (GCMT) (1976 -December 2010), The compiled instrumental earthquake catalog is covered from 1900 to 2013 ( Figure 1). All reported event magnitudes have been converted to moment magnitude by using Scordilis (2006) relationship for mb-Mw and Ms-Mw. Subsequently, all duplicated events have been removed, and earthquake magnitude and location correction have been performed manually to remove any obvious errors. The largest earthquake magnitude 9.2, Great Sumatra-Andaman earthquake in northern Sumatra on 26 December 2004 at 0 7.58 local time.

Magnitude conversion
In the final updated earthquake catalogue, several different magnitude scales are used to define the earthquake magnitude. For example, the 20-s surface-wave magnitude (M S ) and the short-period P-wave magnitude (mb) are commonly used in the data from USGS, ISC, and other international database sources, and the moment magnitude (M W ) is reported in the Global Centroid Moment Tensor catalogue. It is necessary to convert all these different magnitude scales into a single magnitude scale. In this study, the moment magnitude scale is chosen as the single representative scale. Since the accuracy of reported magnitudes is dependent on magnitude definitions, the more reliable magnitude is then preferred for using in magnitude conversion as follows: M w , M S ,mb, and M L Conversions between magnitude scales are made using the equation provided in Table 2 in Ornthammarath et al. (2011). After the magnitude conversion, we merged duplicate entries (from different data sources) into a single entry for each earthquake event.

Declustering
One basic assumption of the adopted seismic hazard assessment methodology is that earthquake occurrences are statistically independent (the Poisson assumption). Therefore, the earthquake catalogue to be used for seismic hazard assessment must be free of dependent events, such as foreshocks and aftershocks. The process to eliminate dependent events from earthquake catalogues is called "declustering". Gardner and Knopoff(1974) declustering algorithm, is chosen for the present study. This approach states that foreshocks and aftershocks are dependent (a non-Poissonian process) on the size of the main event, and these earthquake events need to be removed in accordance with space-and time-windows. Normally, a large main earthquake event leads to larger aftershocks over a larger area and for a longer time. Therefore the time-and distance-window parameters for larger main events are greater than those for smaller events. Declustering eliminates about 60% of the 64,866 events in the catalogue. The final declustered catalogue includes 25,654 earthquake events (4218, 7585, and 13851 events for shallow, intermediate, and deep earthquake, respectively) with M W greater than or equal to 3.0 in the study region from 1900 to 2013.

Catalogue completeness
It is recognized that earthquake data in the catalogue are not complete, and that failure to correct for the data incompleteness may lead to underestimation of the mean rates of earthquake occurrence. The correction can be made by identifying the time period of complete data for prescribed earthquake magnitude ranges. Reliable mean rates of earthquake occurrence for the given magnitude ranges can then be computed from the complete data.
Two methods were employed for completeness analysis of the catalogue: (a) the Visual Cumulative method (CUVI) and (b) Stepp's method. Both algorithms provided a similar result; hence, the former technique was adopted. We divide the study region into three zones, i.e., shallow, intermediate, and deep earthquakes (BG-I, BG-II, BG-III, respectively). The data completeness analysis is carried out separately for each of these zones, and the results are presented in Table 2 in Ornthammarath et al. (2011).

MODELING OF EARTHQUAKE SOURCES
To properly describe the complex earthquake environments in the region, they are modeled as a mixture of background seismicity, subduction area sources, and crustal faults. These are described in more detail below.

Background seismicity model
The background seismicity model represents random earthquakes in the whole study region except the subduction zones. The model accounts for all earthquakes in areas with no mapped seismic faults and for smaller earthquakes in areas with mapped faults. In this approach, it is not necessary to divide the region into many small areas. One large area may be used, but the rate of seismicity is assumed (or allowed) to vary from place-to-place within the area. The rate of seismicity is determined by first overlaying a grid with a given spacing, in the current case 0.10° in latitude and 0.10° in longitude, approximately 10 by 10 km, onto the study region, and counting the number of earthquakes with magnitude greater than a reference value (M ref ) in each grid cell. The rate of seismicity is computed by dividing the number of earthquakes by the time period of earthquake data. The rate is then smoothed spatially by a Gaussian-function moving average and comparing with the observed seismicity. By this approach, the spatially-varied seismicity can be modeled with confidence relating to source uncertainty. In hazard calculations, earthquakes smaller than magnitude 6.0 are characterized as point sources at the centre of each grid cell, whereas earthquakes larger than magnitude 6.0 are assumed to be hypothetical finite vertical or dipping faults centered on the source grid cell. Lengths of finite faults are determined using the Well and Coppersmith (1994) relations. Consecutively, the pre-calculated average source-to-site distance from virtual faults with strike directions uniformly distributed is employed (Petersen et al. 2008).
The whole study region is divided into 6 source zones: BG-I, BG-II, BG-III, SD-A, SD-B, SD-C, (see Fig. 3). The zones SD-A, SD-B, and SD-C are subduction zones, which will be described in detail below. The zone BG-I is a background seismicity zone covering the whole study area excepting subduction zone, and the zone BG-II and BG-III are background seismicity zone for intermediate and deep earthquake except subduction zones. Earthquake data, particularly greater than magnitude 3.0 earthquakes, in BG-I, BG-II, and BG-III could be reliably recorded since 2004 due to the high earthquake detection capability of a fairly dense seismograph network in Malaysia and surrounding region. Hence, the accuracy of the estimated seismicity rate in BG-I can be significantly improved by including small earthquakes (3.0<M W <5.0) in the seismicity rate calculation Furthermore, this activity rate computation is also based on the observation that moderate earthquakes generally occur in areas where there have been a significant number of magnitudes 3 events. On the other hand, in BG-II and BG-III earthquake data with M W > 5.0 can be used for computing the seismicity rate due to the incompleteness of small earthquake data. Nevertheless, a lack of small earthquake data is not a major problem because the seismicity rate in these zones are relatively high; thus, the rate can be reliably estimated from moderate-sized earthquakes. In addition, the overall influence of BG-II and BG-III on the seismic hazard in Andaman sea is lower than that of BG-I. We model the magnitude-dependent characteristic of the seismicity rate in each background seismicity zone by a truncated exponential model (Gutenberg-Richter model): where N(M W ) is the annual occurrence rate of earthquakes with magnitude greater than or equal to M W , and a and b are the Gutenberg-Richter model parameters. The b-value is assumed to be uniform throughout the whole background region. Hence, we used complete earthquake data with magnitude greater than 4.0 in current study area to compute a single regional b-value. The obtained regional b-value is 0.90, and this value is used for BG-I BG-II and BG-III. The a-value varies from place to place within each zone. It is computed by using a grid with spacing of 0.10° in latitude and longitude and is spatially smoothed using a twodimensional Gaussian moving average operator with a correlation distance parameter C (Frankel 1995). Earthquake data with M W > 3.0 and C = 50 km are used for BG-I, while earthquake data with M W > 5.0 and C = 50 km are used for BG-II and BG-III. The correlation distance is chosen based on Frankel (1995) and it is comparable to earthquake location error. Note that at present there are no fixed rules or guidelines to determine an appropriate C value. If the value of C is too small, the resulting smoothed seismicity will be concentrated around the epicenters of past recorded earthquakes.
On the other hand, if the value of C is too large, the resulting smooth seismicity will be blurred and will not reflect the true spatial variation pattern of seismicity. The chosen C values are believed to suitable as the computed smoothed rate 10a values  (presented in Fig. 2) are in agreement to observed spatial pattern of seismicity.
In the truncated Gutenberg-Richter models of BG-I, BG-II, and BG-III, the minimum earthquake magnitude is set equal to 4.5 because earthquakes with smaller magnitude than this are judged not to cause damage to buildings and structures. The maximum (upper bound) magnitude is set to 7.0 for BG-I and 8.0 for BG-II and BG-III to account for the largest earthquake magnitude that have been observed in these zones, as shown in Fig. 2. The average depth used in the model for BG-I, BG-II, and BG-III are 20, 75, and 120 km, respectively.

Subduction zone model
As explained earlier, the megathrust Sunda subduction zone is divided into seven sub-zones based on seismicity characteristics: the Burma zone (SD-A), the Northern Sumatra-Andaman zone (SD-B), and the Southern Sumatra zone (SD-C). Each sub-zone is modelled as a seismic area source with a uniform rate of seismicity (the traditional area source model), and the magnitude-dependent seismicity rate is modeled by a truncated Gutenberg-Richter relation. The geometry and recurrence times of large earthquakes associated with these active tectonic structures are largely based on available paleotsunami and seismic history studies (Jankaew et al., 2011), summarized geodetic data reported in, Berryman et al., (2013).
The calculated Gutenberg-Richter model parameters (a and b values) are shown in Table 3 in Ornthammarath et al. (2011). The minimum earthquake magnitude in the Gutenberg-Richter model is set to 6.5 as the subduction zones are very near to studied area and hence large earthquakes are also important for long period structures in southern Andaman Sea. The maximum magnitude for each zone is set to the maximum observed magnitude plus 0.5 magnitude units. The maximum magnitude for zone SD-B and SD-C is set to 9.2 as the 2004 Sumatra earthquake and the 2005 Nias earthquake

Crustal fault source model
Two different approaches are employed to model the magnitude-dependent characteristic of the seismicity rate of these crustal faults: the Gutenberg-Richter model and Characteristic Earthquake (CE) model. In the Gutenberg-Richter model, a magnitude-frequency distribution for crustal fault model is assumed from the minimum magnitude of 6.5 to the upper-bound magnitude (Mmax). To account for the uncertainty in estimating Mmax, we consider three different cases with Mmax set equal to MC -0.2, MC, and MC + 0.2. The probabilistic weights of 0.2, 0.6, and 0.2 are assigned to these cases, respectively. In each case, theb-value is set equal to the regional b-value of 0.90, and the a-value is determined from the seismic moment rate, which is computed from the fault slip rate.
In the characteristic earthquake model, three characteristic earthquake magnitudes (M C ) are also considered: M C -0.2, M C , and M C + 0.2. The probabilistic weights of 0.2, 0.6, and 0.2 are assigned to these cases, respectively. In each case, the earthquake occurrence rate is computed from the characteristic magnitude and the fault slip rate (to match with the seismic moment rate of the fault). The recurrence interval for the characteristic model is determined from: where is μ shear modulus, 3.0 × 10 11 dyne/cm 2 , L is rupture length, and W is rupture width, u is the fault slip rate, M 0 C is the characteristic earthquake moment, which is calculated from: log(M 0C ) = 1.5M C + 16.05 (3) and the magnitude is assumed to be normally distributed around the characteristic value with a standard deviation of 0.12. The properties and parameters of Sumatra faults considered in this study are shown in Table 2 in Petersen et al. (2007).

GROUND MOTION PREDICTION EQUATIONS (GMPES)
In this study, three Next Generation Attenuation (NGA) models were developed for shallow crustal earthquakes in the Western United States and similar active tectonic regions were applied for background earthquakes in BG-I and for earthquakes from crustal faults in the study region. These NGA models were developed by Boore and Atkinson (2008), Campbell and Bozorgnia (2008), and Chiou and Youngs (2008) during the NGA project. In addition, based on comparison of several different Ground Motion Prediction Equations (GMPEs) with recorded ground motion from interface subduction earthquakes by Ornthammarath et al. (2014), it has been decided that three subduction GMPEs could be used for probabilistic seismic hazard analysis in this region, and probabilistic weights assigned to these models are 0.25, 0.25, and 0.50 for Atkinson and Boore (2003;2008), Youngs et al. (1997), and Zhao et al. (2006), respectively. These weights are relatively consistent with residual of observed recorded data and estimated values of three selected GMPEs. To calculate ground motion for intermediate-and deepdepth earthquakes, the equations of Young et al. (1997) and Atkinson and Boore (2003) are considered with equal weights.

PROBABILISTIC SEISMIC HAZARD ANALYSIS
The PSHA results for southern Andaman Sea, Figure 4 and 5, are presented in terms of seismic hazard maps at 475-and 2475-year return periods at bedrock condition. For southern Andaman Sea, the observed seismicity in and around Sumatra subduction zone and Sumatra faults control the hazard for most considered structural periods. Estimated bedrock PGA near subduction zone at 2475-year return period range between 0.6g to 1.0g; however, for area near coast of Thailand and Myanmar, estimated bedrock PGA at 475-year return period are relatively less intense varied from 0.05g to 0.15g. This is primarily due to its location which is far removed from any major active structure in augmented with low observed seismicity rate of background seismicity. For long structural period ( Figure 5), large part of southern Andaman Sea is subjected by moderated long period ground motion due to lower attenuation rate of long periods. Subduction zone earthquakes contribute high seismic hazard to long period offshore structure in southern Andaman Sea. Estimated bedrock SA (T = 1s) near subduc-tion zone at 2475-year return period range between 0.4g to 1.0g; however, for area near coast of Thailand and Myanmar, estimated bedrock PGA at 2475-year return period are relatively less intense varied from 0.10g to 0.20g. The suitable ground motion records for performing time history analysis of long period structures in southern Andaman Sea should be selected for large earthquake at long distance. In addition, the computed ground motions are comparable to those in Petersen et al. (2007) with minor different in short period ground motion near subduction zones.