Improving urban seismic risk estimates for Bishkek, Kyrgyzstan, through incorporating recently gained geological knowledge of hazards

Many cities are built on or near active faults, which pose seismic hazard and risk to the urban population. This risk is exacerbated by city expansion, which may obscure signs of active faulting. Here, we estimate the risk to Bishkek city, Kyrgyzstan, due to realistic earthquake scenarios based on historic earthquakes in the region and an improved knowledge of the active fault sources. We use previous literature and fault mapping, combined with new high-resolution digital elevation models to identify and characterise faults that pose a risk to Bishkek. We then estimate the hazard (ground shaking), damage to residential buildings and distribution of losses (economical cost and fatalities) using the Global Earthquake Model OpenQuake engine. We model historical events and hypothetical events on a variety of faults that could plausibly host significant earthquakes. This includes proximal, recognised, faults as well as a fault under folding in the north of the city that we identify using satellite DEMs. We find that potential earthquakes on faults nearest to Bishkek—Issyk Ata, Shamsi Tunduk, Chonkurchak and the northern fault—would cause the most damage to the city. An Mw 7.5 earthquake on the Issyk Ata fault could potentially cause 7900 ± 2600 completely damaged buildings, a further 16,400 ± 2000 damaged buildings and 2400 ± 1500 fatalities. It is vital to properly identify, characterise and model active faults near cities to reduce uncertainty as modelling the northern fault as a Mw 6.5 instead of Mw 6.0 would result in 37% more completely damaged buildings and 48% more fatalities.


Introduction
Many highly populous cities across the world are at risk from earthquakes. Whilst earthquakes at plate boundaries are often larger and more frequent, earthquakes in continental interiors can cause more damage and fatalities despite their smaller size (England and Jackson 2011). This can be due to the location of the fault being unknown, yet having a close proximity to urban centres particularly due to recent urban expansion (e.g. the 2003 Bam earthquake, 2010 Haiti earthquake). This risk can be exacerbated by poor building construction or lack of preparedness and resilience in the affected communities, particularly in middle to low income countries, in contrast to more seismically resilient communities like New Zealand or Japan (Crowley and Elliott 2012). In order to assess the potential seismic hazard and risk posed to cities, we need an understanding of the sources of the hazard: locations and geometry of active faults, the earthquakes that could occur on them, the consequent shaking hazard that could result, and the risk to buildings and population that could be exposed to the hazard.
One of the regions of the world that is rapidly deforming with large active faults and consequently has experienced a number of major continental earthquakes is the Tien Shan (Thompson et al. 2002). The Tien Shan is the intra-continental mountain belt of Central Asia (Fig. 1), which is accommodating roughly two thirds of the shortening of the India-Eurasia Collision (20 ±2 mm/yr, Zubovich et al. 2010) (Fig. 1a,b). It has a history of large (> Mw 7) earthquakes (Kalmetyeva et al. 2009), and is a region of low to middle income countries where major cities have expanded significantly since the last major earthquakes have affected the previously smaller settlements (Fig. 1d).
Bishkek, the capital of Kyrgyzstan, sits at the northern extent of the Tien Shan mountains, and is bounded to the south by the major east-west striking Issyk Ata thrust fault (Fig. 1c). The city has been affected by a number of historically documented and instrumentally recorded earthquakes. The most recent notable events were the 1885 Mw 6.5 Belovodskoe earthquake 50 km to the south-west of Bishkek and the much more distant Mw 7.2 1992 Suusamyr earthquake (Fig. 1c), neither of them occurring on the faults that come closest to the city. In just three decades the city has increased in size to a population of 1 million, compared to 660,000 in 1992 (United Nations 2018), and has expanded southwards towards the major Issyk Ata fault. Bishkek is surrounded by potentially damaging faults, and previous seismic hazard assessments have identified high risk (Abdrakhmatov et al. 2003;Bindi et al. 2012;GEM 2018;Silva et al. 2020;Pagani et al. 2020). However, such previous studies on the seismic hazard and risk for Bishkek are limited in number and require an update based upon more recently acquired geological information of active faulting (e.g. Abdrakhmatov et al. 2003) or inclusion of a wider range of scenarios (e.g. Bindi et al. 2012) who consider rupture only of the Issyk-Ata fault).
Bishkek has an increasing population and the city is growing. This means the risk is potentially increasing, due to the growing number of buildings and people, or the construction type of new buildings, or the relative proximity of new developments to the sources of hazard as the city expands along, towards and onto active faults. Previous studies have shown that smaller, proximal earthquakes can cause significantly more damage than larger distal events (Hussain et al. 2020;Amey et al. 2021), which demonstrates the need for detailed fault studies, combined with hazard and risk assessment. b a c d Fig. 1 Setting of Bishkek and the Tien Shan, Central Asia. a GNSS velocities from Zubovich et al. (2010) across the Tien Shan, relative to a fixed stable Eurasia, on a SRTM 30 Plus DEM showing elevation (Becker et al. 2009). b Profile through these GNSS velocities, taking along line X-X' in panel a), with red indicating profile-parallel component of velocity, and blue indicating profile-perpendicular component of velocity c The setting of Bishkek city, with its four districts shown in black and the major identified faults shown in red, notable earthquakes from Kalmetyeva et al. (2009); Chediya et al. (1998) d Bishkek urban extent, with colour indicating build-up in 2014 derived from Landsat images (Pesaresi et al. 2016;Corbane et al. 2018). Over time the city has expanded geographically, in addition to the centre becoming more densely built-up. On a hillshaded and coloured SRTM30 background (Farr et al. 2007) Whilst it is fairly intuitive that distance to a fault rupture is one of the dominant factors in estimates of shaking, the aim of this study is to provide quantitative estimates of the relative potency of different realistic fault rupture scenarios in terms of the potential range of losses and their distribution.
In defining the scenarios, we base these potential earthquake ruptures upon the growing knowledge of the active faulting in the region from recent geological, palaeoseismology and remote sensing investigations (e.g. Landgraf et al. 2016;Mohadjer et al. 2016;Patyniak et al. 2017), and provide evidence for further scenarios from our topography analysis presented here. We additionally provide scenarios based upon past events where our knowledge of them from recent research better constrains their locations and magnitudes (e.g. Ainscoe et al. 2018;Krüger et al. 2018).
These quantitative estimates are achieved through a Scenario-based Seismic Hazard assessment based on a collection of rupture geometries and a set of ground motion models for shallow active crust earthquakes. Additionally, we aim not only to estimate average losses but also to capture the range of potential outcomes based upon uncertainties in the expected ground shaking and vulnerability. The epistemic uncertainty in the selection of the most suitable ground motion model is incorporated through the consideration of various GMPEs, while the aleatory uncertainty in the ground shaking is propagated into the risk results through the sampling of hundreds of ground motion fields, considering the standard deviations of each GMPE. The uncertainty in the vulnerability component is encapsulated through the consideration of a probability distribution for each loss ratio (Sousa et al. 2016). For each scenario, we present not only the average losses, but the standard deviation and range of outcomes to highlight the long tail in the distribution, due to the incorporation of these measures of uncertainty.
In this paper, we use high-resolution satellite-derived Digital Elevations Models (DEMs), including TanDEM-X datasets and a DEM of Bishkek city that we derived from Pleiades tri-stereo optical satellite imagery to identify proximal faults that could potentially fail in earthquakes. We then use the OpenQuake seismic hazard and risk assessment software to compute scenario damage and risk estimates for specific, defined earthquake events Silva et al. 2014). We use OpenQuake and its associated tools due to several reasons including its open-source nature, ability to incorporate the aleatory uncertainties in the ground shaking and vulnerability components, large library of GMPEs, and capacity to estimate both hazard and risk. A major advantage of this approach is that it raises the transparency of the data analysis performed here, increasing the scientific reproducibility, and enabling quicker updates for those that might want to modify or introduce new scenarios for the city. We present the potential shaking hazard (peak ground accelerations) and risk (building damage, fatalities and economic cost) from multiple events, extending the work of Bindi et al. (2011) who calculated risk scenarios on one approximation of the Issyk Ata fault. We have previously highlighted the impact of the choice of fault segmentation and geometry with depth of the estimate of ground shaking losses for nearby Almaty, Kazakhstan (Amey et al. 2021). We test the impact that the choice of the sub-surface geometry, segmentation and position of rupture for the Issyk Ata and Shamsi Tunduk fault immediately south of the city has on our calculation of estimated losses to highlight the potential impact of this epistemic uncertainty.
We assess which districts of Bishkek are most at risk from earthquake shaking, and identify the earthquake scenarios that would cause the most damage if they were to occur. This provides the evidence for which districts of the city are at greatest risk from known faults, and also which parts are more exposed to potentially hidden faults. The aim of such individual scenarios is to present a clear case of outcomes in terms of losses which will aid in the communication of risks to motivate preventative actions in future spatial planning.
2 Description of the methodologies adopted for the identification of ruptures and the calculation of hazard and risk

Identifying fault scarps and building heights from high-resolution DEMs and field GPS measurements
We use DEMs, coupled with previous mapping and literature reviews, to map active faults around Bishkek, in order to identify and parameterise earthquake scenarios for hazard and risk analysis with OpenQuake. We use TanDEM-X DEMs of 12 m resolution (Borla Tridon et al. 2013), as well as a new high-resolution DEMs of Bishkek city created from Pleiades tri-stereo optical imagery using Agisoft Metashape (version 1.6.2). This photogrammetric process identifies tie-points between the pair of images, performs triangulation and produces a dense point cloud of latitude, longitude and elevation, making use of the rational polynomial co-efficients (RPCs) provided with the satellite imagery data. We processed the Pleiades stereo panchromatic images at Ultra High Resolution (i.e. the native resolution of the satellite images) with mild depth filtering. Details of pointcloud availability can be found in the acknowledgements. From the point clouds, we produce DEMs. With these DEMs we take profiles of topography to identify from the surface any evidence of folding to the north east of Bishkek (Fig. 2), over anticlines to the East of Bishkek and within Bishkek itself (Sect. 4.3).
Folding to the north-east of Bishkek has previously been identified (Thompson et al. 2002). In 1910 there was a notable earthquake to the north-west of Bishkek which damaged the village of Georgievka, now Korday/Kordoy in Kazakhstan, on the border with Kyrgyzstan. This event, known as the Georgievskoe earthquake, is reported to be a magnitude 5.6 at 43.00 • N, 74.48 • E and 10 km depth, south-dipping (Kalmetyeva et al. 2009). This is ∼15 km west of the location of the northern fault we model here, indicating there is precedent for earthquakes in this region.
We take profiles over this folding using TanDEM-X topography, as well as GNSS (Global Navigation Satellite System) measurements which we took during a field campaign in 2019. Using a Leica GS18 GNSS attached to a Leica CS20 hand-unit, we mounted the GS18 on a 2 m pole and took measurements every second from a car as we drove along the roads. The GS18 was positioned out the window above the car and accounted for the height in our measurements. We used the Smartlink function, such that the GS18 acted as a rover without a base-station for real-time kinematic operation, with corrections supplied by geostationary augmentation satellites. For the profiles, we take 10 m bins and 2 m swaths on the TanDEM-X along the profile lines shown in Fig. 2a, and 20 m width and 2 m profiles for the SRTM 30 m DEM. The elevation is not easily identified in a hillshaded DEM (Fig. 2a), but in a colourised DEM (Fig. 2b) and in elevation profiles (Fig. 2c), the increase in elevation can be clearly observed. We note the strike of the surface expression of folding may be modified by the Chu river. We ascribe this folding to a fault propagation fold, with uplift caused by a thrusting fault at the base, and map a fault of 10 km length at the base of this increase in elevation, to use as an input source geometry to the OpenQuake Engine for the risk calculations (see Sect. 2.3).
We also took profiles within Bishkek city to identify if faulting had stepped north into the basin, potentially underneath Bishkek city itself, as has been identified to the East in Almaty, Kazakhstan (Grützner et al. 2017;Amey et al. 2021). In order to identify subtle slope changes within a city, it is necessary to use high-resolution DEMs, which here we have created from Pleiades tri-stereo optical satellite imagery at 2 m resolution. However, a DEM of this resolution will include anthropogenic features of the city such as buildings and trees. In order to assess the terrain, we use lastools (Isenburg 2020) to create a digital terrain model (DTM): the ground surface, with features like buildings removed. We classify the highest and lowest points in each 3.5 m cell, and use this to classify the ground and 'non-ground' points, to create the DTM of the ground surface. Once the anthropogenic features have largely been removed, we can then take profiles through the city. We take profiles of 10 m bins and 2 m swaths, and on profiles that show offsets we calculate the gradient of the ground using 1 km length, to remove any noise of inadequately removed buildings in the DTM. The results of these profiles and potential buried faults are discussed in Sect. 4.3.
We use the faults we have identified in satellite DEMs, along with those previously identified and discussed in literature, to define a number of earthquake scenarios to use as inputs for the hazard and risk calculations in the OpenQuake engine.

Exploring earthquake hazard and risk scenario calculations using the OpenQuake engine
The OpenQuake engine is an open-source hazard and risk modelling software, created by the Global Earthquake Model (GEM) foundation Silva et al. 2014). We use OpenQuake's Scenario Calculator (version 3.11.3) to compute the hazard (ground shaking, e.g. peak ground acceleration, PGA), losses (fatalities and economic cost in US dollars) and damage (buildings showing damage) due to specific, defined earthquake scenarios. The scenarios that we use are defined in detail in Sect. 2.3. For each rupture scenario we define the earthquake parameters: fault location, geometry including dip and depths, earthquake magnitude and mechanism. We then use the ground motion prediction equations (GMPEs) of Abrahamson et al. (2014), Campbell and Bozorgnia (2014), and Chiou and Youngs (2014) to calculate the ground shaking as a function of distance from the fault plane. We use these GMPEs as they were formulated for shallowtectonic settings such as this one, and they require distance relative to the closest point on the fault ( R rup ), whether horizontal distance or vertical distance for locations above the fault plane. The GMPEs use estimates of shear wave velocity in the top 30 m ('Vs30') to account for site amplification affects, for which we use the USGS Vs30 estimate using slope as a proxy (Allen and Wald 2009) as we do not have measured velocity values for Bishkek (Fig. 10). To account for uncertainties in the GMPEs, we run each calculation  (Farr et al. 2007). b Zoom in of the profile lines, with colour indicating elevation from the same SRTM30 DEM used in panel a, highlighting the elevation increase to the northeast of Bishkek. Lower panels show the profiles, with orange TanDEM-X (which does not cover the whole area), SRTM30 shown in blue (42 m have been subtracted from the SRTM to match the geoid to ellipsoid), and GNSS measurements for profile D, which follows the orientation of a road we drove along during a field campaign in 2019. Using this elevation increase we have mapped a fault, and where this fault crosses the profile is shown by the red line ▸ 1 3 In order to calculate damage and losses we require exposure information: details of the building types in Bishkek and occupancy. The exposure in Bishkek has been extensively collected through field campaigns and using satellite imagery by Wieland et al. (2012); Pittore and Wieland (2013); Wieland et al. (2014); Pittore et al. (2020). They categorise building attributes including building material and height throughout Bishkek. Figure 3a shows the distribution of exposure points in Bishkek, approximately every 1.5 km. Here, we perform calculations on the residential buildings at night, meaning that all residents are assumed to be in the building at time of the event. Figure 3b shows the distribution of buildings types for the four districts of Bishkek (Leninsk, Pervomaisk, Oktyabrsk and Sverdlovsk).
We include exposure points outside the official district spatial definitions, as the exposure dataset includes many points very close to, yet not within, these districts as currently defined. If these outer exposure points were not included, the exposure dataset would only be a population of ∼580,000 and 14,000 buildings, which is too little for Bishkek: the 2018 UN estimate of population of Bishkek in 2020 is ∼1,000,000 (United Nations 2018). We include outer exposure points in the nearest district using a 5 km buffer, as shown in Fig. 12, which means our exposure dataset for Bishkek has a population of ∼870,000 and 30,000 residential buildings.
For the vulnerability component, we use the fragility and vulnerability functions proposed by Martins and Silva (2020), which used an analytical approach to develop functions that covered the building classes included in the exposure model for Bishkek. This a c b Fig. 3 a Four administrative districts in Bishkek (Leninsk, Pervomaisk, Oktyabrsk and Sverdlovsk) with colour indicating the total number of residential buildings in that district. The white points show the location of exposure points used in this study (Wieland et al. 2012;Pittore and Wieland 2013;Wieland et al. 2014;Pittore et al. 2020). The black lines indicate the outline of the city districts, with the thin lines indicating the expanded extent of each district used in this study in order to fully capture the exposure in the region (see Fig. 12). b Total replacement cost of residential buildings per district, in units of million US dollars. c Total occupants per district. d Residential building type break-down in each district. Each district has unreinforced masonry as the prevalent building type derivation method accounted for a wide range of uncertainties, including the record-torecord and building-to-building variability. The fragility curves (i.e. probability of exceeding a set of different damage states conditional on ground shaking intensity) were converted into a set of vulnerability functions (i.e. probability distribution of loss ratio conditional on ground shaking intensity) using a damage-to-loss model (i.e. relation between the expected loss ratio for each damage state). In this process, the uncertainty in the fragility functions and damage-to-loss models is propagated into the vulnerability functions using the total probability theorem (e.g. Silva 2019). The consideration of these sources of uncertainty is particularly important for the type of analyses presented herein due to the significant heterogeneity of the building stock in Bishkek.

Earthquake scenario selection
We identify a total of nine scenarios: three historical and six hypothetical. This enables us to estimate how much damage these historical earthquakes would cause to Bishkek if they were to happen today, using our exposure model's estimate of the current population. Based on the fault mapping done in this study and those of others (details in following Sections) we identify faults near to Bishkek and estimate consequences to the city if earthquakes were to occur on these today. Figure 4 shows these faults in relation to Bishkek, as well as their surface expression and depth. Each earthquake is described in detail below. All the parameters are given in Table 1. Details on how to access the input files for Open-Quake can be found in the acknowledgements.

Historical scenarios
The following earthquakes are shown with Fig. 4 with purple indicating the estimated length of the projection of the fault to the surface, and the blue dashed lines indicating the assumed outline of the fault plane at depth.

Balasogun Mw 6.4 1475
The Balasogun earthquake is estimated to be a magnitude ∼6.4 (Kalmetyeva et al. 2009), with an approximate epicentre was 80 km south-east of Bishkek. The earthquake

Table 1
OpenQuake parameters used to model eight earthquake scenarios. All scenarios extend from the surface (top depth = 0 km). The 'ID' corresponds to the numbers in 20 Thompson et al. (2002), Mohadjer et al. (2016) destroyed the town of Balasogun and settlements in the surrounding region, including the 11th century Burana fortress (Chediya et al. 1998). We assume this earthquake occurred on the Shamsi Tunduk fault, as no evidence of surface faulting has been found. Kalmetyeva et al. (2009) report an epicentre of 42.6 • E and of 75.2 • N. Given the nature of faulting within the immediate region, in our model we use a purely reverse earthquake event (rake = 90 • ), dipping uniformly at 45 • to the south, extending from the surface to 20 km depth and ascribe it a magnitude of M w 6.4. Belovodsk Mw 6.9 1885 The Belovodsk earthquake occurred on the 2nd of August, 1885 to the south-west of Bishkek (Kalmetyeva et al. 2009;Landgraf et al. 2016). It is believed to have a moment magnitude of ∼6.7 (Bindi et al. 2014) − 6.9 (Kalmetyeva et al. 2009). This is reported to have destroyed the villages Belovodsk and Kara-Balty (Kondorskaya et al. (1982) in Bindi et al. (2014)), causing a total of 54 fatalities and 77 injured. The report of Ignatiev (1886) gives similar values, with 37 fatalities and 43 injured people in Belovodsk, and 17 fatalities and 20 injured people in Karabalty. No fatalities or injuries were recorded in Bishkek at this time.
Using the MI L H (surface wave magnitude) contour lines from Bindi et al. (2014), and the description of the location of the surface rupture in Kalmetyeva et al. (2009), we model this earthquake as rupturing along the mountain range front, on the Issyk Ata fault, with a length of 25 km, with a moment magnitude of 6.9. We use a rake of 90 • , and dip the fault at 45 • to the south from the location and extent of the surface to 20 km depth.

Suusamyr Mw 7.2 1992
The most recent large (> M 7) earthquake to have occurred in Kyrgyzstan is the Suusamyr earthquake on 19th August 1992. This occurred in the Suusamyr region, approximately 100 km south-west of Bishkek, with 50 fatalities, damage and landsliding in the surrounding villages (Ghose et al. 1997), and minor structural damage to Bishkek. The moment magnitude for this earthquake is estimated at magnitude 7.2 Gómez et al. 1997). The surface ruptures have been well mapped (Ainscoe et al. 2018;Grützner et al. 2019) which we use to define the surface rupture.
We model this earthquake as magnitude 7.2 on a 50 km long fault. We model it as dipping at 60 • to the south and a rake of 120 • Gómez et al. 1997), extending from the surface to 20 km depth.

Hypothetical scenarios
Anticline to East 50 km to the east of Bishkek, there are a number of anticlines which may be caused by fault propagation folds. These are potentially along-strike continuations of the Zailisky Alatau fault range to the east around Almaty (Grützner et al. 2017) and are identified in the GEM Global Active Faults Database (Styron and Pagani 2020). We plot profiles over this anticline (Fig. 11) and map a fault based on this scenario. We delineate it as three segments, dipping at 45 • to the south from the surface to 20 km depth. Based on a 70 km fault length we model it as capable of accommodating a M w 7.5 earthquake, with a rake of 90 • .

Chonkurchak fault
The Chonkurchak fault is located south-west of Bishkek. From the west it follows the mountain range front, before continuing east into the mountains, where the Issyk Ata fault follows the mountain range (Fig. 1). We model a scenario in which the 70 km length Chonkurchak fault ruptures as one fault, in a large event. Trenches at Panfilovkoe, located on visible fault scarps likely associated with northward propagation of the Chonkurchak fault by Landgraf et al. (2016), give a paleomagnitude of 6.7 and 7.2/7.3 of historic events, indicating that this fault can fail in large earthquakes that will likely be damaging to Bishkek.
As such, we model the Chonkurchak scenario with an event of magnitude of M w 7.5, dipping at 45 • south, extending from the surface to 20 km depth. Landgraf et al. (2016) found dip-slip motion in their trenching at Panfilovkoe, so we model a rake of 90 • .

Shamsi Tunduk
The Shamsi Tunduk, on which the 1475 Balasogun earthquake may have occurred, is one of the closest faults to Bishkek. It is located south-east of Bishkek and formed the high mountains visible from the city (with the lower hills closest to Bishkek formed by the Issyk Ata fault). We explore the impact of failure of the whole Shamsi Tunduk fault in one event. We model it as a 70 km long rupture extending from the surface to 20 km depth, with a magnitude of M w 7.5, with a dip of 45 • south and a rake of 90 • .

Issyk Ata
The Issyk Ata fault is closest to Bishkek and bounds the city to the south (Abdrakhmatov 1988), separating the Tien Shan mountains to the south from the Kazakh platform in the north. Sections of its 120 km length have been mapped in various studies (Thompson et al. 2002;Landgraf et al. 2016;Patyniak et al. 2017). The Belovodsk earthquake occurred on the western extent of the Issyk Ata fault, and as such may have transferred stress to the rest of the Issyk Ata fault, pushing it closer to failure in an earthquake. As such we model the Issyk Ata fault from the eastern extent of the Belovodsk rupture modelled here (see inset in Fig. 4) and extend it 70 km, as is consistent with an Mw 7.5 earthquake. We model it with a rake of 90 • , dipping at 45 • from the surface to 20 km depth.
Folding to the north Folding to the north of Bishkek has been identified in previous literature (Thompson et al. 2002;Landgraf et al. 2016). In Sect. 2.1 we plotted profiles over this folding (Fig. 2) and suggest it could be caused by a fault-propagation fold. The 1910, Mw 5.6, Georgievskoe earthquake occurred ∼10 km west of this folding (Fig. 1), showing there have been moderate earthquakes in this area. The surface expression of this faulting may have been modified by fluvial erosion of the Chu river, and potentially the fault may strike more East-West in line with the mountain range front than on Fig. 2b, but has been subsequently eroded. In addition to the damage in the building stock, an earthquake here could potentially have other risk implications for Bishkek, as there are several dams located on this elevated ground.
Using our profiles over this fold to constrain its east-west extent, we model a 10 km length fault, with a M w 6.0 rupture occurring on it, dipping at 45 • south-west. We extend this fault from the surface to 10 km depth, with a rake of 90 • . We note that the Chu river flows through here, and it is possible this has removed the trace in some sections. We explore the difference an impact of an alternative fault model with a larger magnitude would make in Sect. 4.5.1.
North-dipping fault to the south Further to the south, the faults in the Suusamyr valley continue further to the east, as identified by (Thompson et al. 2002;Mohadjer et al. 2016). Two prehistoric earthquakes have been identified as ∼ 3 kyr ago, and >8 kyr ago by trenching at locations that did not fail in the 1992 Suusamyr earthquake (Ainscoe et al. 2018), indicating several earthquakes have been caused by south-dipping faults in this valley. There is evidence of a north dipping fault along the northern edge of the valley, as seen in scarps and late Quaternary alluvial and moraine deposits (Fig. 1b Ainscoe et al. 2018), which have not been trenched, for which we use as a basis for this scenario.
To explore how a large event like this would affect Bishkek compared to more proximal faults, we model a fault of 150 km length in three segments, dipping at 45 • north, from the surface to 20 km depth. We model this as a magnitude of M w 7.8 earthquake with a rake of 90 • . An earthquake of this magnitude and distance would be similar to the north-dipping 1911 Chon-Kemin earthquake (Arrowsmith et al. 2017) and its affect on the city of Almaty in Kazakhstan (Amey et al. 2021). Figures 5 and 6 show the hazard and risk results for the different scenarios considered, with totals given in Table 2 and Fig. 7. In Figs. 5 and 6 we have aggregated the damages and losses across the four Bishkek districts (which have been extended, as discussed in 2.2), in order to show the spatial variation across Bishkek. The bar charts show the percentage loss or damage for that district. The damage that we present is the number of completely damaged residential buildings -this measure includes buildings that completely collapse, as well as those that are damaged to such a degree that they must be demolished and rebuilt. The full distribution of building damage states (slight, moderate, extensive, complete) is shown in Figs.13 and 14 and in Table 3. The damage by building type is shown in Table 4, which shows that unreinforced masonry buildings make the up the highest percentage of completely damaged buildings: in our scenarios unreinforced masonry buildings constitute between 80% and 94% of all completely damaged buildings. The replacement cost is the cost of replacing the residential buildings (not contents) in US dollars (2017), and may show a different pattern across the city depending on the average cost of a residential house in that district. The number of fatalities can also show a different spatial distribution, depending on the prevalent residential building in that district (e.g. high-rise flats compared to detached houses). The percentage of fatalities is also significantly lower in each district than the damage, e.g. a district may have 1% of buildings in that district completely damaged but only 0.1% of the people in that district will consequently die. This is because only a small proportion of completely damaged buildings fully collapse, and in addition even those that do collapse are not certain to cause fatalities. Most of the fatality models indicate fatality rates between 10% and 40% in completely damaged buildings, depending on the type of construction in the fragility/vulnerability models used here Martins and Silva (2020).

Illustration of earthquake hazard and risk results
The historical scenarios show lower damage and losses than the hypothetical scenarios. This is due to the larger distance from Bishkek for earthquakes of their size (Mw 6.9 and 50 km away for Belovodsk, Mw 6.4 and ∼30 km away for Balasogun, Mw 7.2 and over 60 km away for Suusamyr). Of the three historical scenarios, Belovodsk would be the most damaging if it were to occur today, with significant damage of 660 ± 700 buildings completely damaged, 330 ± 200 million USD in replacement costs and 130 ± 200 fatalities. In comparison, we estimate 50 ± 100 and 40 ± 100 completely damaged buildings for Balasogun and Suusamyr, respectively, 100 ± 100 million USD and 10 ± 0 fatalities for both Balosogun and Suusamyr. We compare these to observed damages and losses in Sect. 4.1. The uncertainty range is calculated by taking the standard deviation of the 1000 realisations calculated per GMPE.  Fig. 4). Second row shows the Peak Ground Acceleration (PGA), in units of g, which is the mean of the three GMPEs used here (Abrahamson et al. 2014;Campbell and Bozorgnia 2014;Chiou and Youngs 2014). The next three rows show, respectively: damage (number of buildings classified as completely damaged), cost (replacement cost of the residential buildings in US dollars), fatalities, aggregated over the four city districts, which are extended to capture the city exposure better (see Fig. 12). The bar charts surrounding the city districts indicate the percentage damage, cost or fatalities for the four districts

Fig. 6
Results from OpenQuake scenario modelling of the hypothetical earthquake scenarios. Input details are given in Table 1. The top row shows the location of the fault input (for full details see Fig. 4). Second row shows the Peak Ground Acceleration (PGA), in units of g, which is an average of the three GMPEs used here. The next three rows show, respectively: damage (number of buildings classified as completely damaged), cost (replacement cost of the residential buildings in US dollars), fatalities, aggregated over the four city districts, which are extended to capture the city exposure better (see Fig. 12). The bar charts surrounding the city districts indicate the percentage damage, cost or fatalities for the four districts. The scenarios that incur the most damage and losses are the northern fold and Issyk Ata / Shamsi Tunduk scenarios For the hypothetical events, the faults closest to Bishkek cause the largest damage, losses and fatalities. These scenarios are the Mw 7.5 Issyk Ata, the Mw 7.5 Shamsi Tunduk, the Mw 6.0 Northern fold and Chonkurchak faults. We estimate the Issyk Ata event would cause 7900 ± 3500 completely damaged buildings and a further 16,400 ± 2000 slightly, moderately or extensively damaged buildings, 2400 ± 1000 million USD in replacement cost and 2400 ± 1500 fatalities. The northern fold event would cause 2600 ± 1700 completely damaged buildings, 900 ± 500 million USD replacement cost and 600 ± 600 fatalities. The Chonkurchak fault would cause a similar level of damage and losses as the northern fault, despite being further from the city, due to its larger magnitude. We estimate an earthquake on the Chonkurchak fault would cause 2200 ± 1800 completely damaged buildings, 800 ± 500 million USD in replacement costs and 540 ± 600 fatalities.
For the northern fold, this level of damage despite the comparatively small magnitude is due to fault dipping under the city, with a number of exposure points located directly above the fault experience a high degree of shaking in the hanging wall. While the Issyk Ata and Shamsi Tunduk faults do not intersect with the city and the exposure points as much, the larger earthquake magnitude of 7.5 results in peak ground accelerations across the city. This means, the damage, losses and fatalities are anticipated to be spread more evenly across the city districts (Fig. 6). For the northern fault, the majority of damage and losses occur in the north-eastern district (Sverdlovsk) and to a certain extent in the north-western district (Pervomaisk), due to the fault extending under these districts.
We stress that while the Issyk Ata, Shamsi Tunduk, Chonkurchak and northern fault scenarios are the most damaging, even the faults at great distance cause significant damage. The Eastern Anticline and Southern fault are the least damaging hypothetical scenarios explored here, but would still cause hundreds of buildings to completely collapse, resulting in over a hundred fatalities and over $300 million of direct economic losses in replacing and repairing damaged buildings. Therefore, despite their large distance from Bishkek, these are significant events.

Discussion
In our deterministic seismic hazard approach, the aleatory uncertainties of estimating the ground motion fields from specified earthquake ruptures is encapsulated in the implementation of the GMPEs. We additionally have selected three different GMPEs appropriate for shallow-tectonic settings to reduce the impact of relying solely on one subjective choice of model. Additionally, in the risk calculation using vulnerability, we also account for the uncertainty in estimating loss ratios that propagate through to the final results. The distribution of consequence losses calculated then have a very large standard deviation that in many cases is as big as the average loss expected ( Table 2). The compounding of these uncertainties can result in departure from normal distributions with some long tail large outcomes in the estimated losses (Fig. 7). Additional information regarding the impact of these sources of uncertainty in earthquake scenarios can be found in Silva (2016). In this section, we discuss some of the choices in the fault characterisation that would further add to the level of epistemic uncertainty in our results, and highlight where future research is needed to better constrain these parameters that might result in a reduction in the variability of the risk estimates.

Comparison to historical losses
We compare the damage and loss estimates for the two most recent historical earthquakes modelled here to those recorded from these events in the historical records. There is very little known about the Balasogun earthquake in 1475, so unfortunately we have no records to compare our results to. This highlights a limitation with our approach in using the historical catalogue of earthquake ruptures, as records are relatively short for this area compared to other countries such as Italy, in addition to the establishment of permanent settlements being much more recent for this region. The short record compared to the recurrence interval for earthquakes in this region means that the small sample of historic earthquakes will not be representative of the ones that Bishkek will eventually be exposed to. This motivates us to then use the recently gained knowledge of active faulting from geological, palaeoseismic and remote sensing investigations to better sample the possible rupture events in this region that could affect Bishkek.
In the 1885 Belovodsk earthquake, Ignatiev (1886) reported no fatalities or injuries in Bishkek itself, though the villages Belovodsk and Kara-Balty to the west of Bishkek are reported to have been destroyed, with over ∼ 50 fatalities and ∼> 60 injuries (Ignatiev 1886;Kondorskaya et al. 1982), since Bishkek has expanded significantly since 1885 (Fig. 1d shows expansion since only 1975). Because of the increased population and building stock, we estimate the Belovodsk event would cause 660 ± 700 buildings completely damaged and 130 ± 200 fatalities if it were to occur today. This estimate is for Bishkek alone, and as in the 1885 event it is likely that the villages further west, closer to the rupture, will be worse affected.
The 1992 Suusamyr caused damage and fatalities in the surrounding villages (Ghose et al. 1997) but only minor damage in Bishkek. Our models estimate 30 ± 100 completely damaged buildings and 10 ± 0 fatalities. These figures are low due to the distance of Bishkek from the rupture (>60 km), but Bishkek could potentially see slightly more damage and loss than incurred in 1992 due to the increase in urban population as shown in Fig. 1d. Along-strike eastward propagation of rupture due to increased stresses from Suusamyr is explored in one of the hypothetical scenarios.

Comparison to other literature
The seismic risk for Bishkek from the Issyk Ata fault has been estimated by Bindi et al. (2011). Unlike in this study, they model the Issyk Ata fault as striking due east, dipping at 50 • south. The western extent of their modelled fault is 74 • longitude, whereas our modelled Issyk Ata fault's western extent is 74.25 • , as we assume the segment that failed in the Belovodsk event is unlikely to rupture again (inset in Fig. 4). Bindi et al. (2011) estimate the total number of collapsed buildings as 22,200 (>3× higher than estimated here) with a further 44,410 damaged, and the fatalities are 16,600 (>7× higher than estimated here). We attribute this difference in part to Bindi et al. (2011) using an exposure model with >77,000 buildings, compared to the ∼30,000 we use here, though their exposure model has a similar population -850,000 in total, and we assume a population 870,000. We also attribute this difference to an over-estimation of the fragility and vulnerability by Bindi et al. (2011), as they find over 91% of buildings would experience some kind of damage. This is a higher percentage of damaged buildings than has been observed in similar events, for example Gorkha in which the magnitude is smaller and the vulnerability may potentially be higher.
With Bishkek identified as at high seismic risk, there have been a number of studies investigating the feasibility of early warning systems in Bishkek (Picozzi et al. 2013;Bindi et al. 2015;Parolai et al. 2017). For the Belovodsk earthquake, Picozzi et al. (2013) estimate that Bishkek would have between 5 and 9 s lead-time (time available in which to take protective measures). The closest scenarios to Bishkek and amongst the most damaging, the Issyk Ata and northern fault, would have less than 5 s lead time, with an assumed seismic network of 19 strong-motion stations (Parolai et al. 2017). In contrast, an earthquake on the southern fault described in this study would have 5-10 s lead time Parolai et al. (2017). This lack of preparation time for proximal earthquakes will likely further exacerbate the damage caused, as even with significant improvements to early warning systems the lead time is still very small.

Potential faults within Bishkek city
Cities built along mountain range fronts are at risk from faults stepping into the basin, for example as seen in: Landgraf et al. (2016) in which fault scarps are likely associated with northward propagation of the Chonkurchak fault west of Bishkek; in Tehran (Talebian 2004); or as in seen in Almaty city in Kazakhstan (Grützner et al. 2017;Amey et al. 2021).
In order to investigate whether there are faults in the city, we use the Pleiades tristereo DEM, as discussed in Sect. 2.1. Whilst there are a number of locations on separate profiles that show an offset or change, there are three locations that consistently show changes across multiple profiles (Fig. 8). In some profiles this may be: an offset, a change c and d show the profile W and C, respectively, with the measured offsets visible, but in the opposite sense to that shown by the Issyk Ata fault in panel b. Note vertical exaggeration in panels b, c and d in gradient, a localised area of uplift, or there may be no change. There is also one at ∼ 8 km, which we attribute to the railroad. Figure 8 shows the location of the profiles we have taken, and shows offsets for lines W and C with Table 5 showing the status of the gradient for each profile, either recording an offset, or identifying a change in gradient, localised elevation increase, or no change. However, the offsets that we observe do not show the sense of motion one would expect for northward propagation of the Issyk Ata fault. If the Issyk Ata fault had propagated into the city, we would expect to see southward dipping faults, with uplift on the southern side (foot wall) and subsidence on the northern side (hanging wall). But the profiles in Fig. 8 show the opposite: uplift on the northern side. As such, it is not clear if these locations represent faults propagating into the city. These could potentially represent the onset of folding. If these coherent offsets in the alluvial fan were back-thrusts, we would expect to see frontal thrusts or folding, which we do not observe. The city grid is also orientated N-S and E-W, and fault features that we are looking for would be orientated E-W, making it harder to distinguish tectonic features from urban developments.
Additionally, not all the profiles clearly show offsets or even a change in gradient. This may be due to the DEM and its processing (e.g. removing buildings to create a digital terrain model), due to the city and anthropogenic change, or the possibility of variable uplift along-strike of a potential fault.
We recommend further work in the form of GNSS measurements at these locations and field observations, in order to ascertain the cause of these offsets, and whether they could represent faulting in Bishkek city. The scenarios we have discussed here show that earthquakes with comparatively small magnitudes (e.g. the northern fault modelled at Mw 6.0) can cause significant damage and loss, if the faults are close to Bishkek. As such it is important to rule out whether or not there is faulting within the city.

District and building classes most at risk
In Figs. 5 and 6 we plot the percentage damage, cost and fatalities in each district. Which district is most affected is highly correlated with the location of the rupture; for proximal events, the district closest to the rupture shows the highest magnitude and percentage damage and losses, for example in the Belovodsk, Northern fold, Issyk Ata and Shamsi Tunduk scenarios. But for distal earthquakes, at which distance the peak ground accelerations experienced across the city are fairly uniform at the scale plotted in 5 and 6, the two northern districts (Sverdlovsk in north-east, Pervomaisk in north-west) show consistently higher percentage damage than the southern districts. This is seen in the Balasogun, Suusamyr, Eastern Anticline, and Southern fault scenarios. This is similar to previous studies: Picozzi et al. (2013) estimate a fairly uniform distribution of damage index across the city for a Belovodsk scenario, but the highest level of damage is expected only in the north of Bishkek. Bindi et al. (2012) also ascribe higher intensities in the northern part of Bishkek than in the south.
In contrast to what might be expected given these findings, Fig. 3d shows that Sverdlovsk and Pervomaisk have a higher percentage of reinforced masonry, and a lower percentage of unreinforced masonry than the other two districts, meaning that we might expect these districts to perform better under similar ground shaking as reinforced masonry is less likely to be damaged than unreinforced masonry. This suggests the cause of this slight increased percentage damage in the northern districts could be due to local effects or building type. The southern part of town is built on alluvial fans, whilst the north is built on increasingly finer fluvial sediments, which is reflected in the Vs30 estimates (Fig. 10), and likely contributes to the percentage damage spatial pattern seen here. In addition, the Sverdlovsk district has a high percentage of people living in high-rise flats, as is reflected in Fig. 3, as it has a similar number of residential buildings to Leninsk district in the south-west: 8600 residential buildings for 280,000 occupants is ∼ 33 people per building in Sverdlovsk, compared to 9,000 residential buildings for 215,000 occupants is ∼ 24 people per building in Leninsk. We use a topographic slope proxy approach to estimate the average shear wave velocity to 30 m as we do not have Vs30 measurements for the city of Bishkek to better characterise the site conditions. Whilst this approach has been found to act as a reliable proxy in other areas when compared against in situ measurements (Allen and Wald 2009), this adds an additional level of uncertainty to our loss calculations that could be improved through individual site investigations.

Northern fold
Here, we have modelled the northern fault as a magnitude 6.0 earthquake on a fault of 10 km length, that extends to 10 km depth. But from the profiles shown in Fig. 2 and observations of the folding, it would be reasonable to extend this fault to 18 km length, and consequently extending to 15 km at depth, further under the city assuming an equi-dimensional rupture. A fault of this size would equate to a larger magnitude 6.5 event (Wells and Coppersmith 1994), and we model this rupture to see the impact on damage and losses.
We find that increasing the length, depth and magnitude of this northern fault would result in 37% more completely damaged buildings, 31% more economic cost and 48% more fatalities. The reason for this significant increase is due to both the fault extending further under the city, and the larger magnitude resulting in higher magnitude of ground shaking across the city. The characterisation of this scenario has the largest epistemic uncertainty of all those considered here as we have weak constraints on the likely fault geometry. Here, the estimated losses are not the largest of those explored but changes to our assumptions would have the greatest impact on this model. This shows the importance of properly characterising faults close to cities, as the surface expressions of faults can help inform assumptions about the fault geometry at depth and maximum magnitude of earthquake that could occur on it.
In both these scenarios we have assumed this fault dips at 45 • from the surface to depth. If the fault were instead to dip and shallow to decollement beneath Bishkek, the damage and losses could be significantly increased, as has been shown in Almaty scenario risk models (Amey et al. 2021). As the estimated losses are most sensitive to the depth to faulting beneath the city, this is one of the key parameters to constrain in future assessments of faulting in the area in order to reduce the uncertainties in our modelling.
We also note that this fault is directly underneath a number of dams in the region; two large (several kilometres across) and a series of smaller ones (< 1 km). Further work is needed to investigate the impact of earthquakes on these dams, including whether supply of water delivery would be disrupted due to pipe damage.

Issyk Ata and Shamsi Tunduk faults
The Issyk Ata is the closest fault to Bishkek, and forms the boundary of the southern expansion of the city, with many people living within a few kilometres of this fault. Here, we have modelled an earthquake on the Shamsi Tunduk fault and separately an event on Fig. 9 Comparison of earthquake scenarios on the Shamsi Tunduk (f), Issyk Ata (k) and a scenario in which the Issyk Ata fault fails as part of the Shamsi Tunduk rupture (a). The second row shows peak ground accelerations experienced across Bishkek (districts shown in black outline, fault shown in red at surface, blue at depth). Row three shows the number of completely damaged residential buildings aggregated for each district. Row four shows the economic replacement cost aggregated for each district, in millions of US dollars (2017). The bottom row shows the number of fatalities aggregated for each district the Issyk Ata fault, but it is possible that an event on the Shamsi Tunduk fault would also trigger failure on the eastern segment of Issyk Ata fault.
We model an earthquake scenario on the Issyk Ata and Shamsi Tunduk faults, in which the Issyk Ata fault dips at 45 • from the surface down to 3 km, then an approximately flat (5 • ) ramp extends from 3 km depth to 3.5 km depth, before then extending down to 20 km depth dipping at 45 • (indicated by dashed lines in Fig. 9a). This is consistent with cross sections by Bullen et al. (2001); Thompson et al. (2002), that suggest the Issyk Ata only extends to shallow depths in Neogene sediments before joining the Shamsi Tunduk fault.
Our results show that the Issyk Ata fault failing as part of a rupture on the Shamsi Tunduk fault would cause more damage and losses to Bishkek than rupture purely on the Shamsi Tunduk (Fig. 9), as the Issyk Ata is several kilometres closer to Bishkek. Our estimates suggest that the Issyk Ata also failing in an event of the same magnitude on the Shamsi Tunduk (Mw 7.5) would incur 39% more completely damaged buildings, 36% more economic replacement cost and 50% more fatalities. This is more similar to the loss and damage seen in our Issyk Ata scenario. This variation in losses (despite all being modelled as having the same magnitude) demonstrates the sensitivity of the result to the choice of the fault structure selected. This emphasises the need for good structural cross sections to constrain the segmentation and geometry of faulting, in addition to identifying which splays are active in faulting through detailed palaeoseismic investigations.

Conclusions
In this study, we have used satellite DEMs to characterise and identify faults that could pose hazard and risk to Bishkek. We have used the open-source software OpenQuake Engine to calculate damage and losses to Bishkek for specific, defined earthquake scenarios. We present the total completely damaged buildings, economic replacement cost and fatalities for Bishkek as a whole and its individual districts, as well as percentages for each district.
Bishkek is surrounded by active faults, and this study has shown earthquakes on many of these would potentially be very damaging for Bishkek. Of the scenarios modelled, we find that the Mw 7.5 Issyk Ata, Mw 7.5 Shamsi Tunduk and Mw 6.0 Northern events would cause the most damage due to their proximity to Bishkek, despite the northern fold being the smallest earthquake modelled. Our calculations estimate the Issyk Ata event would cause 7900 ± 3500 completely damaged buildings, a further 16,400 ± 2000 slightly, moderately or extensively damaged buildings, 2400 ± 1000 million USD in replacement cost and 2400 ± 1500 fatalities. The northern fold event would cause 2600 ± 1700 completely damaged buildings, a further 12,500 ± 1,700 damaged buildings, 900 ± 500 million USD replacement cost and 600 ± 600 fatalities. We recommend future work building on these results could incorporate these scenarios and the relative faults into a hazard model that would allow the calculation of hazard and risk also considering the occurrence frequency of each event.
Our results highlight the need to properly characterise the location of faults close to major cities, including understanding their geometry and interaction with other faults, Fig. 10 Vs30 estimates used in this study, using USGS's topographic estimation (Allen and Wald 2009). Black outline shows the Bishkek districts, and white dots show the exposure points at which we have estimates of building type and number. In the GMPE calculations, the nearest Vs30 estimate is used to calculate ground motions at each exposure point particularly in Bishkek where the Issyk Ata fault and fault-propagation fault under the northern folding. We also recommend further work be done, including GNSS field campaigns within Bishkek city, in order to rule out active faulting within Bishkek city. As the northern fault scenario has shown, even a small earthquake occurring close to Bishkek can be extremely damaging.
The Issyk Ata and Shamsi Tunduk faults to the south and folding to the north dominate Bishkek's landscape; they are a constant, stark reminder of the earthquake risk to Bishkek. Understanding these faults, their geometry at depth and how they interact, as well as how they may connect if at all to faults north of Bishkek, is vital to understanding the seismic hazard and risk to Bishkek.

Fig. 11
Top -The Eastern Anticline fault modelled in this paper shown in red, with the outline of the Bishkek districts in blue. Black lines A-E show the location of the profiles. Bottom -Topographic profiles from SRTM30 and TanDEM-X, in the locations shown in the top panel. Red gives an indication of where the fault modelled in this paper crosses the profile, though note the angle of this fault is not to scale Fig. 12 In this paper, we have expanded the districts to include exposure points that don't fall within the city district outlines, but which contribute significantly to the exposure. Panel a shows the Bishkek district outlines (black), and exposure points included in this (green) and those excluded (orange). A number of those excluded are very close to the district outlines, and the estimated 580,000 is too low for Bishkek. Panel b shows the original city districts (thick black) and our expanded city districts (thin black) in which we have padded the districts by 5 km, and again showing those exposure points included (green) and excluded (orange). The occupancy estimate for this is 870,000, which is more reasonable to population estimates (United Nations 2018). Panel c shows which exposure point is included in which district  Table 3 Estimated residential building damage from the earthquake scenarios modelled here. In the main text we show completely damaged residential buildings, i.e. those damaged to such a degree that they must be demolished and rebuilt, but as is shown here there would be a significant number of slightly, moderately or extensively damaged buildings. These all contribute to replacement cost Damage  Table 4 The break-down of residential building types that were completely damaged in each earthquake scenario. For each scenario, the first column represents the total of each building type that showed completely damage, and the second column indicates the percentage of the total. Unreinforced masonry significantly contributes to the total of completely damaged buildings, ranging from 80 to 94% of all completely damaged buildings Building class