Seismic hazard assessment of Western Coastal Province of Saudi Arabia: deterministic approach

Seismic hazard assessment is carried out by utilizing deterministic approach to evaluate the maximum expected earthquake ground motions along the Western Coastal Province of Saudi Arabia. The analysis is accomplished by incorporating seismotectonic source model, determination of earthquake magnitude (Mmax), set of appropriate ground motion predictive equations (GMPE), and logic tree sequence. The logic tree sequence is built up to assign weight to ground motion scaling relationships. Contour maps of ground acceleration are generated at different spectral periods. These maps show that the largest ground motion values are emerged in northern and southern regions of the western coastal province in Saudi Arabia in comparison with the central region.


Introduction
A major reason for anxiety arises from the increase of seismic vulnerability of most urban structures especially in developing countries (Re 2000). The most important element in the seismic design of the structure and its vulnerability assessment against the earthquake damage is ground motion vibrations. These parameters are evaluated through probabilistic seismic hazard assessment (PSHA) or deterministic seismic hazard assessment (DSHA) approaches (Waseem et al. 2013). The DSHA depends on selected earthquake scenario, specific ground motion probability, and closest distance to the site of interest. However, PSHA deals with numerous scenarios for earthquake and ground motion probabilities (Lin and Baker 2011).
In DSHA investigation, a maximum credible earthquake is supposed to occur at the closest distance to the site of interest. Meanwhile, the occurrence of earthquake probability during the specific period of exposure is neglected (Tavakoli et al. 2013). It must be noted that there is no commonly accepted deterministic seismic hazard analysis approach that is applicable to all parts of the world as it is experienced in assorted ways in different parts of the world. The traditional DSHA methodology delineates the seismic source or sources that might affect the site of concern and then computes the maximum possible earthquake magnitude for each of these sources. By considering each of these maximum earthquakes placed at a location that put the earthquake at the minimum possible distance to the site, the ground motion is expected, mostly, applying an empirical attenuation relation.
Previous studies, about seismic hazard assessment along Western Coastal Province of Saudi Arabia are limited and focus only on PHSA [e.g., Al-Amri (2013) and Al-Arifi et al. (2013)]. Recently, there are few published researches related to DSHA that are applied to some local cities located within Western Coastal Province of Saudi Arabia [e.g. Almadani et al. (2015)]. Yet there is no publication which accounts for the mapping of DSHA that could cover whole Western Coastal regions of Saudi Arabia. The deterministic approach is subject to the epistemic and aleatory uncertainties (Stepp et al. 2001). In this research, the epistemic uncertainties are treated by taking alternatives for the ground motion prediction equations, which in concern link several different assessments of the ground motion. The aleatory uncertainties are related to various input parameters used to describe the seismicity and the ground motion prediction equation.
In the present study, DSHA is carried out along the whole Western Coastal Province of Saudi Arabia (Fig. 1). The grid spacing for regional studies is usually 50 km 9 50 km or half degree (latitude or longitude) e.g., (Deif et al. 2013;Mohamed et al. 2012). A grid of 50 by 50 km (about 155 points) is selected along the eastern coast of Red Sea. For each point, the DSHA is evaluated with associated uncertainty and adopted by a variable logic tree (essentially due to change in distance between points of the grid). In current study, the used methodology is proposed by Campbell (2005) and applied by Deif et al. (2013Deif et al. ( , 2009. This methodology deals with the uncertainty problem and provides ground motion at different percentile levels at various spectral periods. This provides great flexibility to the engineer to select the appropriate input ground motions for selected site. 2 Crustal structure and regional seismicity Lithospheric structure and composition of the Arabian Plate are constructed using a broad choice of geophysical data (Stern and Johnson 2010). The Arabian Plate has a varying crustal thicknesses range. It varies from 22 to 53 km across the plate from west to east (Fnais et al. 2013). The variation is at modest amount in central Arabia ranging from 32 to 46 km in west to 35-50 km in east (Al-Damegh et al. 2005;Rodgers et al. 1999). The crustal thickness approaches 43 km at east of the Central Arabian Magnetic Anomaly (CAMA) and 38 km thick in western part of CAMA. There is a considerable crustal discontinuity in lateral direction beneath the CAMA, causing velocity increment around 0.2 km/s in northeast (Gettings et al. 1986).
The dominant mafic composition of the Arabian Shield is the reason for high crustal velocities in lower crust. The studies of entire shield by surface wave designated the Moho's depth to be at 41-46 km (Mokhtar 2004;Mokhtar et al. 2001). The crust composition, which is inferred from velocity structure, indicates lower mafic crust overlying by upper felsic crust, having low velocities (Stern and Johnson 2008).
The Arabian plate seismicity is clustered and controlled along its major tectonic borders, Red Sea rifting, Dead Sea transform fault, Zagros fold and thrust belt, Biltis thrust, and Makran subduction zone. Seismic activities in the vicinity of the study area are mainly controlled and narrowed along the Red Sea rifting, Gulf of Aqba, and Gulf of Aden (Fig. 2). The seismic activity is concentrated along Red Sea axial trough (Al-Malki and Al-Amri 2013).

Methodology
The estimation of ground motion parameters is carried out through seismic hazard assessment (SHA) at a site of interest for seismic design (Hashemi et al. 2013). SHA analysis depends on seismic activity and attenuation relations (Dowrick 2009). However, DSHA is restricted to a specific earthquake scenario.
The methodology for the current study followed Campbell (2005) procedure for the DSHA which is described in context to the study area as below: 1. Identification of the seismogenic sources having possible impact on Western Coastal Provinces of Saudi Arabia. 2. Calculation of maximum earthquake (M max ) for each seismic source zone using earthquake catalogue and Gutenberg and Richter relationship. 3. Selection of a specific set of earthquake magnitude and distance scenarios. 4. Selection of appropriate set of attenuation laws (ground motion scaling relationships). 5. Incorporation of uncertainties using logic tree. 6. Calculation of the median ground motion for each scenario. 7. Delineate the largest median value among calculated median ground motions. 8. Calculate fractile (percentile) largest median value.
In the current study, the logic tree is utilized to capture uncertainty. Thus the following structure of current methodology for DSHA is build up; identification of seismic sources, calculation of maximum credible earthquake magnitude M max , selection of appropriate ground motion attenuation equations, and finally building logic tree to assign weight to GMPE.

Identification of seismogenic sources
The seismic source delineation is a major key parameter in hazard assessment (Vipin and Sitharam 2013). The seismogenic source characterizations follow historical and recent seismicity, seismicity pattern, seismogenic potential of active faults (Meletti et al. 2008). The seismogenic source model is identified and delineated by consideration from ancient times and current seismicity, tectonic and geological setting, crustal tomographic studies, crustal heat flow measurements, and current hazard studies (Al-Amri 2013; Al-Arifi et al. 2013;Al-Damegh et al. 2005; Al-Malki and Al-Amri 2013; Ares 2010; Burkhard and Grünthal 2009;Pailoplee et al. 2010). In current study, seismogenic source model by Rehman (2016) is utilized to carry out hazard assessment process (Fig. 3).

Earthquake catalogue
One of the most important products of seismology is earthquake catalogue which provides a broad dataset in earthquake events. This can be used in various analyses associated with seismicity and seismotectonic, hazard assessment, and physics of earthquake. The hazard parameters are determined well if the catalogue has longer time coverage (Woessner and Wiemer 2005). A catalogue is a basic requirement for seismicity analysis in space-time volume, and in SHA (Gupta et al. 2012;Leonard et al. 2011).
In current study, the earthquake catalogue is compiled for spatial region extending from 15°to 35°N and 29°to 47°E and include events from magnitude 3 and above from 827 to 2013 AD. This catalogue is used to characterize seismicity of study area. It is compiled by combining information from different sources which include  Poirier and Tahir (1980), Badawy (1999), Ambraseys et al. (2005a, b), Ambraseys et al. (1995Ambraseys et al. ( , 2009), Guidoboni et al. (1994), and Guidoboni and Comastri (2005). The completeness magnitude (M c ) is defined theoretically as the lowest magnitude at which 100 % of earthquakes are detected in space-time volume. The precise estimation of M c is critical because higher values of M c lead to undersampling, and too low values are erroneous. Mainly catalogue-based and network-based techniques are applied for M c estimation (Mignan and Woessner 2012). Magnitude of completeness is a basic requirement to model seismicity in an area. The maximum curvature technique is mainly used (Wiemer and Wyss 2000) for the completeness magnitude.
In current study, M c for each seismogenic source zone is calculated using maximum curvature technique. The Fig. 4 is an example of M c calculation for one zone.

Gutenberg and Richter relationship (G-R)
The seismicity of the seismic source zone is described by the means of famous recurrence relationship termed as G-R relationship.
Log N ¼ a À bM; N represents the earthquakes of specific magnitudes (M) or larger per year, a is activity rate and defines the intercept of the above equation at M = 0. The factor b is the slope which depicts the comparative proportion of small and large magnitudes.
The Gutenberg and Richter (1944) relationship introduces an impractical supposition in which the largest size possible earthquake in any zone being studied, is unrestrained and unconnected toward seismotectonic setting. Kijko andSellevoll (1989, 1992) extended the Gutenberg-Richter equation from data that contain large historical events and recent observation with that of different quality and heterogeneity. In current study, the b values are calculated using Kijko and Sellevoll (1992) approach for each zone (Table 1).

Maximum magnitude (M max )
Another most important parameter beside recurrence parameters is maximum expected magnitude (M max ). The SHA is strongly influenced by the choice of M max . The M max can be estimated using either deterministic or probabilistic approach. Deterministic approach comprises empirical regression relationships between the magnitude and various tectonic and fault rupture parameters, however the probabilistic approach involves extreme value statistics (Gupta 2002). In current study, M max for each seismogenic source zone is calculated using methodology proposed by Kijko (2004) based on observed seismicity. This method is designed for the calculation of the maximum earthquake magnitude, M max , for a seismogenic source zone. This technique is based upon broad equation for the estimation of M max , which is very flexible. This equation is applicable to the following three different scenarios: • Distribution of earthquake magnitude follows doubly truncated Gutenberg-Richter relation. • Moderate deviations occur from the Gutenberg-Richter relation for empirical magnitude distribution. • There is no specific form of magnitude distribution assumed.
The first two solutions of the generic equation are parametric and the third solution non-parametric. The M max values calculated for each seismogenic source zone are described in Table 1.

Ground motion scaling relationships
The ground motion scaling relationships are functions of seismological parameters which include earthquake magnitude, source to site distance, and site conditions (Atkinson and Boore 1995; Peruš and Fajfar 2009;Yazdani and Kowsari 2013). Joyner and Boore (1981) postulated functional form to derive these empirical relationships. The geometrical spreading for all distances is accounted in this functional form (Ambraseys et al. 2005a;Ambraseys et al. 1996;Boore et al. 1997;Sabetta and Pugliese 1987;Smith et al. 1996). The selection of GMPE is accomplished preliminarily from a list of available equations.
The selection and rejection criteria for a GMPE in this study are based upon Cotton et al. (2006) guidelines. The key consideration for GMPE selection/rejection is as described below: • Irrelevant tectonic regime, • Insufficient dataset, • Model is not yet published in peer-reviewed journal, • Unsuitable frequency range for engineering applications, and • Inappropriately regression analysis or regression coefficients misjudged.
Six different ground motion scaling relations are selected in accordance to Cotton et al. (2006) guidelines, which fulfill the tectonic conditions of the study area. These models include Abrahamson and Silva (1997), Campbell and Bozorgnia (2003), Sadigh et al. (1997), Atkinson and Boore (1995), (Boore et al. 1997) and one next generation attenuation equation by Campbell and Bozorgnia (2008).

Logic tree
The logic tree comprised several branches which portray confidence on different input parameters. Weight is assigned to each parameter according to the confidence level on each input (Abrahamson and Bommer 2005;Bommer 2002;Bommer and Scherbaum 2008). In current study, logic tree framework is utilized to incorporate several GMPE. The GMPE are assigned different weights in accordance with applicability of each model with respect to source-site distance. The source to site distance is divided into three ranges 0-200, 200-600, and 600-1000 km. Maximum weight is assigned to next generation attenuation model (Fig. 5).

Results and discussion
The DSHA is carried out with EzFrisk 7.52 TM program and ground motion acceleration values are calculated in cm/s 2 at bedrock. The deterministic response spectrum is generated for MEAN and 0.84 fractile. The maps for MEAN and 0.84 fractile at peak ground acceleration (PGA), 0.1, 1, and 2 s spectral periods are generated (Figs. 6,7,8,9,10,11,12,13).
The ground acceleration values are extremely high in northern part of the study area. The possible explanation for such acceleration is because of close proximity of northern part of study area to seismogenic sources 19, 21 22, and 23 (Fig. 3). The ground acceleration decreases toward central part where moderate acceleration values are present. However, the moderate values can be categorized into two distinct behaviors for accelerations (Figs. 5,6,7,8). The values for ground acceleration are comparably low for shoulder area along the Red Sea coast as compared to its eastern part. The controlling source for this area is Zone 7 and Zone 15, which have low seismicity. However, the eastern parts of central area are controlled by higher seismicity zones than western part (Zones 6, 10, and 14). At southern part of the study area high ground acceleration behavior is observed. The controlling sources for this part of the study area are Zones 1, 2, 4, and 5.

Conclusions
The deterministic approach actually considers the worst case ground motion. The worst case ground motion affects the design and cost of building. These ground motion spectra generated are used as essential input parameters for the upgradation of design parameters for local regulatory  requirements in accordance with international standards. The design parameters are dependent on the local geology and local site soil conditions. In the current study we calculated ground motion at bedrock. However, the effect of soil must be taken into consideration for local site-based design parameters. The DSHA technique is utilized to determine maximum expected ground acceleration for western coastal regions of Saudi Arabia. This research concludes that northern and southern part of western Saudi Arabia are more prone to higher seismic risks compared to central coastal areas. A detailed seismic risk assessment based on site response analysis is recommended for whole coastal regions. The top priority shall be given to northern and southern coastal regimes. This detailed risk assessment will facilitate qualitative decision making for redundant industrial systems and post-earthquake recovery plans.