Probabilistic seismic hazard assessment of Bishkek , Kyrgyzstan , considering empirically estimated site effects

It is well known that variability in the surface geology potentially leads to the modification of earthquake-induced ground motion over short distances. Although this effect is of major importance when seismic hazard is assessed at the urban level, it is very often not appropriately accounted for. In this paper, we present a first attempt at taking into account the influence of the shallow geological structure on the seismic hazard assessment for Bishkek, Kyrgyzstan, using a proxy (Vs30) that has been estimated from in situ seismic noise array analyses, and considering response spectral ratios calculated by analyzing a series of earthquake recordings of a temporary seismic network. To highlight the spatial variability of the observed ground motion, the obtained results are compared with those estimated assuming a homogeneous Vs30 value over the whole urban area. The seismic hazard is evaluated in terms of peak ground acceleration (PGA) and spectral acceleration (SA) at different periods (frequencies). The presented results consider the values obtained for a 10% probability of exceedance in 50 years. The largest SA estimated considering a rock site classification of the area (0.43 g) is observed for a period of 0.1 s (10 Hz), while the maximum PGA reaches 0.21 g. When site effects are included through the Vs30 proxy in the seismic hazard calculation, the largest SA, 0.67 g, is obtained for a period of 0.3 s (about 3.3 Hz). In terms of PGA, in this case the largest estimated value reaches 0.31 g in the northern part of the town. When the variability of ground motion is accounted for through response spectrum ratios, the largest SA reaches a value as high as 1.39 g at a period of 0.5 s. In general, considering site effects in the seismic hazard assessment of Bishkek leads to an increase of seismic hazard in the north of the city, which is thus identified as the most hazardous part within the study area and which is more far away from the faults.


Introduction
Probabilistic seismic hazard analysis (PSHA), first formulated by Esteva [1967], Cornell [1968] and Merz and Cornell [1973], is a methodology that allows one to estimate the exceedance likelihood of various levels of ground motion caused by earthquakes at a given location within a given time period.In particular, the goal of PSHA is to quantify the probability of exceeding a certain ground motion level at a site given all possible earthquakes that can affect it within a reasonable distance range.Since Cornell [1968], several modifications have been proposed to account for specific effects that can affect the final assessment.Amongst these, particular attention has been dedicated towards the importance of the local variability of ground motion due to the shallow geological structure.Campbell [1997] proposed to consider the depth of the basement rock when modeling long-period site response, whereas Boore et al. [1997], among others, suggested the use of Vs30 (the average S-wave velocity over the upper-most 30 meters, based on travel time from the surface to a depth of 30 m) as a proxy to account for site effects in seismic hazard assessment.More recently, Campbell and Bozorgnia [2008] and Boore et al. [2008] proposed to estimate 1D site effects based on Vs30 and PGA values through relationships that account for non-linear soil behavior.Campbell and Bozorgnia [2008] additionally attempted to account for 2D-3D effects by modeling the basin response.Note that the above mentioned studies only broadly average the site effects in to the ground motion prediction relationships.Petersen et al. [1997] accounted for site effects in the PSHA of three southern California counties by using ground motion prediction equations representative of three generic site conditions (alluvium, soft-rock and rock-site) that were valid at the regional level, as deduced from geological maps.Steidl [2000] suggested adopting corrections to rocksite-ground motion prediction equations to be used for PSHA at the regional level by estimating regional amplification factors determined by averaging ratios between observed and predicted PGAs and 5% damped SA at 0.3, 1.0 and 3.0 s periods.These amplification factors were therefore dependent on the site-class and the amplitude of ground motion.
Only a few studies have focused on improving probabilistic seismic hazard assessment at the local level.Bazzurro and Cornell [2004] suggested integrating the 1D site response estimated by numerical simulations, while also accounting for uncertainties in the subsoil structure, directly into the PSHA.Erdik et al. [2005] evaluated the seismic hazard for Bishkek accounting for a first-order approximation of the potential site amplification effects, quantifying the building vulnerability and, finally, evaluating the urban earthquake risk.Parolai et al. [2007] used numerical simulations to estimate the site effects in terms of response spectrum ratio for the Cologne (Germany) area.These site effects, along with their uncertainties, were included in hazard analyses by modifying the ground motion prediction equations with respect to a rock site.
In this paper, we carry out a PSHA in Bishkek, the capital of Kyrgyzstan.Bishkek is located in the middle of one of the largest tectonic depressions of the northern Tien-Shan, the Chu basin, over thick Quaternary and Tertiary deposits.It follows that any seismic hazard assessment for the town, which could be used as the basis for a rigorous earthquake risk assessment, must take into account the effect of the shallow geological structure on the spatial variability of earthquake ground motion.Here, we will examine the spatial ground motion variability by considering both the shear wave velocity structure below the city, and the response spectra ratio derived by analyzing the earthquake recordings collected by a urban seismological network installed in Bishkek for a few months [Parolai et al. 2010, Ullah et al. 2013].We use in the first case a standard approach based on the use of Vs30 as a proxy for site effects, while in the second, we follow an approach similar to that of Parolai et al. [2007].Due to the particular sedimentary cover existing in Bishkek (thick sediments up to several kilometers with high shear-wave velocities, especially in the south, of more than 600 m/s) the first approach is expected to be more appropriate for a high-frequency-related parameter such as the PGA.The approach based on the calculation of the response spectra ratio is expected to provide a more accurate description of the ground motion variability for the long to intermediate spectral periods, given that the shorter ones are affected by site effects at the reference rock-station.The obtained results will be discussed in term of the spatial variability of seismic hazard within the urban area, with the limitations and advantages of the followed approaches examined.

Geological and seismic tectonic description of the area
The Chu Basin within which Bishkek is located is bounded by the Kyrgyz range to the south and the Chu-Lli Mountains in the north [Bullen et al. 2001].Below the area of Bishkek, the Paleozoic basement depth is expected to generally increase from the north (~1 km) to the south (~3 km).Seven thick Tertiary formations are overlying the Paleozoic bedrock.Quaternary sediments with a thickness of 200-300 m overlie the Cenozoic deposits.It is worth mentioning that to the south of the study area, there are outcrops of alluvial material made up of rubble and gravel with a thickness of between 15-40 m.The same material can be found within the urban area along the major rivers.In the southern part of Bishkek, alluvial gravels, rubble and sandy material, with a thickness ranging between 25 and 50 m, constitute the shallow geological layers.Lower Quaternary sediments made up of rubbly-beach gravel with claysand lens and break stone beach gravel outcrop in the northern part of the city.In the northernmost part of Bishkek, loess sediments with a thickness of 30-40 m in the east and 60-70 m in the west constitute the shallow geological layers [Parolai et al. 2014].

Earthquake catalog processing
In this study, the earthquake catalog was selected starting from the homogenized catalog developed within the EMCA project for the whole Central Asia region [Mikhailova et al. 2015].The whole catalog, which includes mainly instrumental and historical events starting from around 1800 (only 80 events dating from before 1800 are included) was examined, selecting only "crustal events" (hypocentral depths less than 50 km) and magnitudes (rounded to the 0.1 decimal place) greater than 4.0.Therefore, M 4.0 serves as a minimum magnitude for the following calculations.The catalog was de-clustered using the AFTERAN program [Musson 1999] and a Gardner and Knopoff [1974] distance window.The algorithm implemented in the AFTERAN program is a modification of the Gardner and Knopoff [1974] algorithm which uses a moving time window rather than a fixed one.In the case at hand, a moving time window of 20 days was adopted.More information on the Central Asia catalog processing can be found in Ullah et al. [2015], this issue).
In order to carry out the seismic hazard assessment for Bishkek, two superzones (defined on the basis of seismo-tectonic and geological data) have been selected from the seismic zonation model developed within EMCA [Ullah et al. 2015].Figure 1 shows the distribution of the events located within these superzones after their selection from the de-clustered Central Asia catalog.This new catalog is rich in moderate to large magnitude events, and in particular it includes 12 events with magnitudes greater than 7.While most of the seismicity is concentrated in the south of Kyrgyzstan, in the Tien Shan range, the larger magnitude events (Mw>7) are located to the north of the analyzed area.They includes the MLH 8.3, 1889 Chilik earthquake, which is one of the largest known historic intraplate reverse faulting events, and the MLH 8.2 1911 Chon-Kemin earthquake [Bindi et al. 2013[Bindi et al. , 2014]].The latter is the strongest known event in the Tien-Shan region for which instrumental recordings are also available.The moment magnitude evaluated for the Chon-Kemin earthquake using the available recordings is Mw 7.7 [Storchak et al. 2013, Mikhailova et al. 2015].
The completeness of the selected catalog was assessed by the Stepp method [Stepp 1973], which involved dividing the catalog into bins of 0.5 magnitude width, starting from magnitude 4.0 [Ullah et al. 2015].The analyses are carried out considering time intervals of 5 years, and used the Hazard Modelers' Toolkit of OpenQuake.
The Gutenberg-Richter (GR) parameters i.e., the a-and b-values, are estimated from the de-clustered catalog after completeness analyses using the Weichert [1980] method.This method estimates the maximum likelihood of b-values for grouped magnitudes bins and allows for unequal periods of observation.
Finally, it is worth noting that the catalog is dominated by earthquakes with reverse faulting mechanisms (as indicated by the global CMT catalog; Ullah et al. [2015]), while a smaller number of events appear to display strike slip and normal faulting mechanisms.In particular, in the superzone 7 (Figure 2), 90% of the events show a reverse faulting mechanism and only 10% are strike slip.In superzone 8 (Figure 2), 70% of the earthquakes are reverse fault, 20 % strike slip and 10 % normal faulting events.

Source model
The area sources around Bishkek (Figure 2) were selected from the regional area source model for Central Asia used in EMCA [Ullah et al. 2015].The source zones (being a subdivision of the superzones mentioned above) have been defined considering the distribution of seismicity, the location and strike of the known fault systems existing in the area and a supervised (expert judgment) approach.The sources were selected in order to cover an area within a radius of 200 km from Bishkek.
The parameter related to the completeness of the catalog was derived (as mentioned above) from the superzones and, they were then attributed to the related area sources.The GR parameters a and b for each zone were estimated independently.However, if the number of events are too small (less than 20) or the range of magnitude is too narrow, the estimated a, and b values are not considered to be reliable and the values calculated while analyzing the data for the related superzone are instead used.In this case, the b value of the respec-tive superzone is used, while the a value is calculated using the Weichert method [1980].The results of the analysis, shown in Figures 3 and 4, indicate similar completeness periods (for the same magnitude range) for the two analyzed superzones.We also note from The maximum magnitude was assigned to each zone based on the value estimated for the superzone they belong to.In this case, magnitude values of 7.5 and 8.3 were considered for zones 7 and 8 (Figures 1 and 2), based on the largest known earthquake that occurred in the area, plus the addition of 0.5 based on expert's opinion.The completeness and GR parameters for each of the considered area sources are summarized in Tables 1 and 2, respectively.Parolai et al. [2010] estimated the spatial variability of ground motion in Bishkek using standard spectral ratios (SSR) and the horizontal to vertical ratio (HV) of  the Fourier spectra of earthquake recordings collected at 19 seismological stations installed within the city.Furthermore, Parolai et al. [2010] estimated the horizontal to vertical spectra ratio of seismic noise (NHV) recorded at nearly 200 sites within the urban area of Bishkek.The results allowed a map of the variation in the fundamental resonance frequency within the study area to be developed.Ullah et al. [2013], taking advantage of these results, applied clustering and correlation analyses to the SSR and the NHVs to improve the spatial resolution of site effects in Bishkek.In particular, they identified areas sharing similar site responses in terms of Fourier spectra ratio.Finally, Parolai et al. [2010] and Pilz et al. [2013] estimated the S-wave velocity structure below different sites in Bishkek by means of the analysis of seismic noise collected by micro arrays of seismological stations.They also showed that, using standard relations to obtain Vs30 from the slope of topography [Wald andAllen 2007, Allen andWald 2009] provides misleading values for site response for the Bishkek area.

Site effects
Using the data mentioned above, the influence of the site response in determining the hazard level was studied considering both a proxy for site effects (Vs30) and by modifying the ground motion prediction equations by considering the standard response spectra ratio of empirical data.

Vs30 approach
The ground motion prediction equations chosen for assessing the PSHA in Central Asia within EMCA [Ullah et al. 2015] can account for specific Vs30 values.The Vs30 values used in this study were obtained from in-situ measurements of S-wave velocity derived by the analysis of microarray of noise data [Parolai et al. 2005, Boxberger et al. 2011].Array measurements were carried out ad-hoc, considering the results of Ullah et al. [2013], at sites representative of each area sharing similar site response, following the results of cluster and correlation analysis.
Figure 5 shows the location of the array sites together with a simplified geological map of the area.In general, Vs30 decreases from south to north, consistent with a change in the characteristics of the shallow geological material [Parolai et al. 2014].In the northern part of Bishkek, where the Quaternary sediments mainly consist of silt or loess, the shear wave velocity (Vs30) is low, with values around200 m/s.In the central transition zone, where the sedimentary cover is mainly composed of pebbles and gravels, the Vs30 is of the order of 650 m/s.while in the southern part of the study area, where the shallow surface material consist of boulders, pebbles, and gravel, values as high as ~850 m/s are found.Note that at the station used as the reference for the following site response analysis, Vs30 was estimated to be about 900 m/s.

Response spectrum ratios
As an alternative approach, the site response to be incorporated in the PSHA was also estimated in terms of response spectrum ratios with respect to the reference site [Parolai et al. 2007[Parolai et al. , 2009]].First, the response spectra (5% damping) were calculated for the earthquake recordings of each station of the temporary seismological network.Before calculating the response spectra, the earthquake recordings were de-trended and deconvolved for instrumental response.Second, their ratios with respect to the response spectrum at the reference station were computed.The average of all the response spectra ratios of the stations in each cluster is then found.In agreement with Parolai et al. [2010], the station BI04, located in the Kyrgyz range south of Bishkek, was chosen as a reference site.However, as noted in Parolai et al. [2010], this station, although being the only one located on an apparently rock site, is affected by amplification of ground motion at frequencies higher than 2-3 Hz (0.3 s -0.5 s) and with a maximum at nearly 5 Hz.For this reason, the following analysis was carried out only considering the amplitude of the response spectra ratio (as coefficients to be used in the seismic hazard assessment) over the frequency range 0.2 Hz to 2 Hz (0.5 s to 5 s).Figures 6, 7 and 8 show the average response spectra ratios for all stations.These figures highlight the spatial variability of the response spectra ratios over the study area.For example, in Figure 6, the stations located in the north of Bishkek, where the average Vs30 is as low as 213 m/s, show the largest amplifications.In general, the highest values of amplification are occurring over a broad frequency band.Figure 7 depicts the results obtained for the stations located in the middle of city where the average Vs30 is 650 m/s.These stations show lower average amplifications, mainly affecting frequencies lower than 0.5 Hz and higher than 1 Hz. Figure 8 shows the results obtained for the stations installed in the southern-most and south-east sectors of the city, where the average Vs30 is 862 m/s.These response spectra ratios show a comparatively lower amplification over a broad frequency band.

PSHA including site effects
In this study, site effects are accounted for in the PSHA by considering directly in the ground motion prediction equation (GMPE) the Vs30 values estimated at each site.This leads to the GMPE, which are valid for a rock site, to include the response spectra amplification factors [Bazzurro et al. 2004, Parolai et al. 2007].The Boore and Atkinson [2008] GMPE were used in this study which is also recommended by GEM in active shallow crust.The Boore and Atkinson [2008] GMPE is given in the form: where the terms F M , F D , and F S represent the magnitude scaling, the distance function, and site amplification, respectively.Y is the ground motion parameter, R JB is Joyner-Boore distance metric, M is the magnitude, d T is the period-dependent standard deviation and f is the fractional number of standard deviations.In particular, the site coefficient is represented by considering a linear and non-linear amplification term: The non-linear term F NL depends on the input peak PGA.The linear term depends on the V S30 through the  (3) is described as: where ln(T i ) is the period-dependent average response spectrum ratio for each site.In order to scale the coefficients to a common reference site (consistent with that used for the response spectra ratio calculation) in Equation ( 4), the V S30 used in the linear site specific term is fixed to 900 m/s, independent of the location.Consistent with the PSHA carried out at the regional level [Ullah et al. 2015], the ground motion aleatory variability is truncated at 3 sigma.The PSHA for the urban area of Bishkek was carried out in terms of PGA and spectral acceleration (SA) for different periods.

Results
The PSHA assessment for the area of Bishkek is now shown for the cases when considering an hypothetic homogeneous rock site condition, and for when site effects are introduced.
Figure 9 shows the maximum and minimum limits of the SA with 10% probability of exceedance in 50 years when considering uniform rock site conditions for Bishkek.The values are very close and there is only a slight increase of the ground motion from north to south (Figure 10), very likely related to the dominance of the shaking produced by the larger number of earthquakes sources in the south.The maximum hazard observed for this model is represented by a SA of 0.43 g at 0.1 s (10 Hz) and a PGA reaching 0.21 g in the southern outskirts of the city.
Figure 11 shows the same maximum and minimum hazard as Figure 9 with a 10% probability of exceedance in 50 years estimated using the model including Vs30.The maximum hazard values are 0.67 g for the SA at 0.3 s and 0.31 g for the PGA Note that different from what was obtained when rock-site conditions are considered, the largest values are found in the northernmost area of the city (Figure 12).This is the area where, due to the low S-wave velocities, the highest amplification of ground motion was estimated (see Figure 6).
Similar trends in the results are observed when site effects are accounted for through the response spectra  ratio coefficients.That is, the minimum and maximum SA curves show (Figure 13) larger values than those calculated for a rock site condition (Figure 9) and the area affected by the largest hazard becomes the northernmost one (Figure 14).However, in this case, the maximum hazard calculated for SA at 0.5 s reaches values up to 1.39 g, much larger than that obtained by considering Vs30 only.Note that, in accordance with what was stated above, the calculation for shorter period, and therefore for the PGA, was not carried out.

Discussion and conclusions
In this paper, starting from the results of previous site effects studies, we carried out a PSHA analysis for Bishkek (Kyrgyzstan) that accounts for the modification of ground motion over short distances.Both the Vs30 values used in the GMPE and the response spectra amplification factors introduced into the GMPE were derived from empirical data.
The major effect in including the site response (or a proxy for it) into the PSHA is the re-prioritization of the urban territory in terms of hazard level.In fact, when a generic rock site condition is used, the most hazardous areas are the southernmost ones with a PSHA OF BISHKEK, KYRGYZSTAN  slight decrease of hazard to the north.This result reflects the distribution of the seismicity and the majority of the seismogenic structures with respect to the town.On the contrary, introducing into the PSHA model parameters that describe the effect of the soil response not only masks this effect, but completely inverts the hazard trend in the urban area.
However, the level of hazard estimated by the two proposed procedures disagrees in the absolute values, especially for periods less than 1 second.This might be due to the following reasons: (1) Vs30 is a parameter representative mainly of the very shallow structure that is expected to influence especially the high-frequency part of the ground motion (and therefore the PGA).It might therefore not be the most appropriated parameter in a large basin where the site response at intermediate and long periods might be affected by the shape of the basin itself and where the site response can affect a wide frequency band [Pilz et al. 2011].The ground motion could therefore be underestimated for the medium to long periods (although unlikely to be affected by strong nonlinear behavior).
(2) On the other hand, the response spectral ratio is estimated using weak motion data.While the approach that was followed here can be appropriate for areas with moderate seismicity, it might lead to an overestimation  of the ground motion since it is not accounting for the soil's nonlinear behavior.While this might be a relatively minor problem in the southern parts of the city where the soil mainly consists of boulders and pebbles grading slowly to gravel toward the center of the urban area, it can be a problem for the northernmost part of Bishkek where silty deposit and loess prevail.
It follows that, although this preliminary study showed the importance of accounting for soil response in the quantitative PSHA of Bishkek, further studies are necessary in order to calibrate the procedure for including realistic site effects.The acquisition of in-situ geotechnical parameters, the analysis of borehole data and the combination of empirical data with site responses estimated from non-linear analysis will allow us to improve the seismic hazard and risk estimates for Bishkek, while taking advantage of the track indicated by this study.

Figure 1 .
Figure 1.De-clustered homogenized catalog for Central Asia.Different color and sizes represents different magnitudes ranges as shown in the legend [Ullah et al. 2015].

Figure 2 .
Figure 2. Area source models for Bishkek.The areas sources are from the EMCA seismic zonation model, from which the superzones have been defined.The numbers in red color indicate the numbering of superzones, while the black color indicates the numbering of area sources.
Figures 3 and 4 that the obtained a and b values are very similar, being close to 4 and 0.8, respectively.

Figure 5 .Figure 6 .
Figure 5. Array locations (black stars) with Vs30 values (blue color).The different colors of the stations (single station noise measurements and permanent stations) represent different clusters in terms of the similarity of site effects from SSR (modified from Ullah et al. [2013])

Figure 7 .
Figure 7. Same as Figure 6 except for locations in the middle of the city (see Figure 5).

Figure 8 .
Figure 8. Same as Figure6, except for the southern and south-east sectors (see Figure5).

Figure 9 .
Figure 9.Estimated seismic hazard variation within the city for rock-site conditions (Vs30 = 900 m/s), for 10% probability of exceedance in 50 years.

Figure 10 .
Figure 10.Seismic hazard in terms of 10% probability of exceedance in 50 years for rock-site conditions.The hazard is shown for PGA and spectral acceleration at different selected period in terms of g.For a detailed description of Geological map of the city in the background, see Figure 5.

Figure 11 .
Figure 11.Estimated seismic hazard variations within the city considering site effects based onVs30 variation, for a 10% probability of exceedance in 50 years.

Figure 12 .
Figure 12.Estimated hazard for 10% probability of exceedance in 50 years considering the distribution of Vs30 values.The level of hazard is shown for PGA and spectral acceleration at different selected periods in terms of g.

Figure 13 .
Figure 13.Estimated seismic hazard variations within the city considering site effects estimated from response spectrum ratios, for 10% probability of exceedance in 50 years.

Table 2 .
Summary of the area sources parameters for the study region.*area sources for which the number of events are less than 20.