Probabilistic seismic risk assessment framework: case study Adapazari, Turkey

While earthquakes can have a devastating impact on the economic growth and social welfare of earthquake prone regions, probabilistic seismic risk assessment can be employed to assess and mitigate such risks from future destructive events. In a previous study (Sianko et al. in Bull Earthq Eng 18:2523–2555, 2020), a probabilistic seismic hazard analysis (PSHA) tool based on the Monte-Carlo approach, was developed to predict the seismic hazard for high seismicity areas. In this study, a seismic risk assessment framework is developed incorporating the previously developed PSHA tool, with vulnerability functions based on various damage criteria, exposures and casualty models. Epistemic uncertainty is addressed using logic trees and distribution functions. The developed seismic risk assessment framework can estimate human and economic losses for particular return periods using an event-based stochastic procedure. The framework is applied to a case study area, the city of Adapazari in Turkey. Seismic risk assessment is carried out for different return periods to identify the most vulnerable areas of the city. The verification of the developed seismic risk framework is performed by comparing the predicted seismic losses to those observed during the 1999 Kocaeli earthquake that severely affected the city of Adapazari. The results of the study indicate that while overall predictions for extensive and complete damage states demonstrate strong correlation with the observed data, accurate risk predictions at the district level are not achievable without microzonation studies.


Introduction
Earthquakes can have a devastating impact on economic welfare and resilience of communities, particularly in developing countries. Destructive social and economic consequences of earthquakes on structures and society have been witnessed following several major events (such as Turkey (1999), Haiti (2010), Tohoku (2011) andNepal (2015) earthquakes). The rapid urbanisation of earthquake prone areas in the last few decades makes seismic risk assessment essential in understanding the likely exposure and mitigating its effects. However, the development of seismic risk models for loss estimation is a challenging process due to the numerous parameters involved in the process and their uncertainties.
Commercially available risk assessment tools are normally used by insurance and reinsurance industries. Often these tools are presented as "black boxes" and the user interference is limited to the pre-defined procedures and input parameters (Bommer et al. 2006). As a result, region specific modifications to the hazard and vulnerability models are difficult to implement. Moreover, assumptions and uncertainties adopted in these commercial tools cannot be controlled by the user (Bommer et al. 2006) and the methods used for the conversion of input to output are not generally transparent (Musson and Winter 2012). In recent years, there have been substantial initiatives to develop open source catastrophe modelling platforms with transparent data inputting and outputting procedures and open data standards (e.g. OASIS loss modelling platform (https:// oasis lmf. org/), Aon Impact Forecasting (https:// www. aon. com/)). Moreover, probabilistic seismic hazard models that are available to the public at the regional or country level have been developed to use in seismic risk calculations (e.g. SHARE project (Woessner 2015) and Unified Hazard Tool by USGS (https:// earth quake. usgs. gov/ hazar ds/ inter active)). For example, SHARE model can be used to calculate hazard for required return periods using OpenQuake engine (https:// platf orm. openq uake. org). On the other hand, the hazard results from these models are, in general, available for a limited number of regions or countries for which they were developed and cannot be easily integrated in risk calculations. In addition, these models only provide data for specific soil conditions. Recently, an open access uniform European Seismic Risk Model (ESRM20) has been developed as part of Horizon 2020 EU SERA project (Silva et al. 2020). Open-source earthquake risk software OpenQuake (https:// platf orm. openq uake. org), which also uses event-based probabilistic seismic hazard and risk calculations similar to this paper, is used in ESRM20. However, state-of-art seismic risk studies are required to verify and calibrate the seismic risk models developed within ESRM20. The estimation of damage or losses from past events can provide an opportunity to compare the estimated and observed risks, hence providing valuable cross-check for ESRM20 developers.
Seismic risk assessment requires reliable data on (i) the seismicity, geology and tectonic settings of an area of interest to quantify seismic hazard, (ii) the location and distribution of the building stock and estimates of number of their inhabitants to develop building and population exposure models, and (iii) the characteristics of structural systems of the buildings and their expected performances to estimate structural vulnerabilities/fragilities. Data required for the estimation of the impact of earthquakes on structures and their inhabitants at large scale (e.g. city level or country level) may not be available and/or easily accessible particularly in underdeveloped or developing countries where urbanisation rate is high. To address this issue, global vulnerability assessment procedures using available data can be used. Riedel et al. (2014) developed a data mining method based on Association Rule Learning (ARL) to correlate basic building characteristics (such as number of stories, building age) with the European Macroseismic Scale EMS98 vulnerability classes of structures (Grunthal, 1998) for a given earthquake intensity. As an extension of this work, ARL and support vector machine (SVM) methods were complemented with remote sensing data (e.g. satellite images and aerial photographs) to produce a vulnerability maps for Grenoble and Nice (France) (Riedel et al. 2015). ARL method was also used to perform seismic vulnerability assessment of Constantine, Algeria (Guettiche et al. 2017a). The proposed method gave better prediction of vulnerability when construction material was added to the building features considered in the analysis (e.g. in addition to the number of stories and age of buildings). In a following study, vulnerability proxies created for Constantine, Algeria, (Guettiche et al. 2017a) were applied directly to historical centre of Skikda, Algeria, using available building attributes (e.g. construction period, number of floors, roof shape and material of buildings) (Soltane 2022). In the same study, spatial distribution of expected economic and human losses in Skikda were also predicted based on the models developed by Guettiche et al. (2017b) for the city of Constantine, Algeria using earthquake data from Mediterranean region. An integrated approach, which combines data mining (SVM and ARL), remote sensing, GIS-based mapping, and vulnerability index methods, was developed by Liu et al. (2019) to perform macroseismic vulnerability assessment of the city of Urumqi, China. Although the methods based on data mining are very easy to apply at urban scale as they don't require detailed building information, these are very crude methods and provide rough estimates of relative distribution of damage. Better prediction of seismic risk can be made by using detailed data on the building stock. The study conducted by Konukcu et al. (2017), in which remote sensing techniques (e.g. areal and satellite images) were used to collect data on the number and age of buildings in Istanbul in GIS format, can be given as example to show how to collect detailed data on buildings.
This work aims to develop a practical yet comprehensive probabilistic seismic risk assessment framework (a new computational platform) using readily accessible data on hazard, vulnerability and exposure. The proposed framework is based on stochastic Monte-Carlo (MC) procedure, that generates synthetic earthquake catalogues by randomizing key input parameters (such as earthquake locations, fault length, maximum earthquake magnitudes, a and b parameters of the Gutenberg-Richter relationship, and ground motion prediction equations). The main advantage of the procedure is that uncertainties in input parameters can be addressed with distribution functions in an efficient manner (Musson 2000). Moreover, logic trees with weightings for each branch can be easily employed within the procedure if required. The main drawback of the procedure is that it can become computationally expensive with both increasing complexity of the model and desired level of accuracy. In the developed framework, fault source zones and background seismicity are considered in the seismic hazard model. Appropriate fragility functions based on various damage criteria and ground motion intensity measures are selected and converted to vulnerability functions using a consequence model to find mean damage ratios (MDRs). While the exposure model is generally obtained from the detailed census data, a practical procedure is proposed to collect building stock data by mapping building footprints from satellite images and gathering data from remote street view surveys when census data are not available.
To demonstrate the capability of the developed seismic risk assessment framework, the city of Adapazari in Marmara region (Turkey) is selected as a case study area. The city is located in a high seismicity area and was previously hit hard by earthquake events in the 20th century. This work contributes to the development of country-specific disaster risk profiles which helps to achieve targets identified within the Sendai Framework for Disaster Risk Reduction (https:// www. undrr. org/ publi cation/ sendai-frame work-disas ter-risk-reduc tion-2015-2030). The results of the study are presented in the form of seismic loss curves and seismic risk maps for Adapazari. Scenario earthquake similar to the 1999 Kocaeli earthquake is employed to compare damage predicted using the developed seismic risk tool in this work with observed damage after that event.

Probabilistic seismic risk procedure
There are numerous seismic risk assessment studies performed for different regions of the world using available engines or frameworks (e.g. Chaulagain et al. 2015;Silva et al. 2015). In this work, a new seismic risk framework is proposed with a general procedure for calculating mean damage ratio using MC simulations. A brief discussion on seismic hazard, fragility functions, consequence and exposure models, which are components of seismic risk calculations, are presented in this section.

Seismic hazard model
Assessment of seismic hazard is an essential component of seismic risk analysis. In a previous study (Sianko et al. 2020), a MC based PSHA was developed as a part of the seismic risk assessment framework proposed in this work. MC simulations are used to generate synthetic earthquake catalogues to represent future seismicity. One of the main advantages of MC procedure over conventional PSHA is its capability of considering the spatial variability with intra-event residuals and efficient way of treating aleatory uncertainties. Conventional PSHA pioneered by Cornell (1968) lacks this advantage and as a result underestimates the total loss value at high return periods (Jayaram and Baker 2009). On the other hand, the spatial correlation of the intra event residuals was not considered in this work as it increases computational complexity. Also, Crowley et al. (2008 and reported that the spatial correlation of the intra event residuals does not influence the average annual loss and does not have a significant impact on the losses for large-scale risk assessment. The detailed information on the MC based seismic hazard procedure employed in this work can be found in Sianko et al. (2020).

Fragility/vulnerability models
Fragility curves can be used to predict the probability of exceedance of certain limit/damage states for a given intensity measure value. It is common for engineers to use simple damage states such as: slight, moderate, extensive and complete damage (e.g. HAZUS-MH MR1, 2003). Consequence models that relate the cost of loss to the rebuilding cost for a given damage state, can be used to convert a set of fragility curves into vulnerability curves to predict economic losses (Erdik 2017;Kohrangi et al. 2021a). Vulnerability models are generally defined as MDR conditioned on ground motion intensity level. MDR can be calculated by integrating consequence models with the proportions of the building stock corresponding to each damage state.
MDR represents the ratio of the cost of repairing a structure, or group of structures, to its replacement cost. It can be used to estimate loss by multiplying it with the economic value of the structure. In the proposed framework, MDR is determined for individual events generated from synthetic catalogues by calculating the ground motion intensity at the site and an appropriate vulnerability model. The worst-case scenario from all seismic sources affecting the site of interest is considered as the annual maximum outcome. This step is applied for all simulations and the results are merged in a single list. The probability of exceedance of certain MDR value can be found by sorting all annual outcomes in descending order and by finding the Nth MDR value in the sorted list. For the desired return period, N can be calculated using the following equation: The sorted list of the obtained MDR values can be plotted against its annual frequency of exceedance to result in a loss curve. Figure 1 shows the flowchart of the procedure for finding probability of exceedance of MDR using MC simulations procedure.

Exposure model
Exposure models include useful information on the location, number of occupants and replacement costs of buildings, in addition to their vulnerability class. Up-to-date exposure data is often unavailable due to the rapidly altering built environment, including aging, particularly in developing countries. Normally, the exposure models rely heavily on the national housing census. Censuses are normally taken every 10 years and are performed at administrative division resolution. The efficiency of the data collected by each country is not consistent, which makes the development of a global exposure model a challenging aspect of risk analysis (Silva et al. 2018).
The main goal of the exposure model is to obtain a layer of uniform resolution across the study area with spatially distributed structures that are classified according to the selected building taxonomy. The taxonomy for the characterization of the exposed building stock and the description of its damage should be compatible with the fragility/vulnerability relationships that will be considered in the risk assessment process (Erdik 2017). For estimating economic and social losses, an exposure model might need to contain additional data about the estimated replacement cost of the structures and expected number of occupants depending on the time of a day. Rapid survey procedures can be performed by utilisation of satellite imagery and from volunteered data in cases where the census data is outdated or missing vital data required for the building exposure model (Wieland et al. 2015).
In addition to housing censuses, there are publicly available sources of data about housing in different countries, e.g. United Nations Housing Statistics (https:// unsta ts. un. org/ unsd/ demog raphic-social/ sconc erns/ housi ng) and the World Housing Encyclopedia (http:// www. world-housi ng. net). Moreover, there are ongoing projects such as Prompt Assessment of Global Earthquakes for Response (PAGER, http:// earth quake. usgs. gov/ earth quakes/ pager) and Global Earthquake Model (GEM, www. globa lquak emodel. org) that have an objective to develop global building inventory databases.

Casualty assessment
Following earthquakes, high casualties are normally observed in highly populated urban areas, where the density of destruction is the highest (Tong et al. 2012). This is particularly important as the number of fatalities is strongly dependent on the number of buildings that collapse or are extensively damaged (Coburn et al. 1992;Feng et al. 2013;So and Spence 2013). There are two main approaches for estimating casualties after an earthquake event.
(1) N = 1 Return period × Catalogue length × Number of simulations + 1 1 3 Fig. 1 The proposed MC-based earthquake risk assessment procedure to calculate mean damage ratio (MDR)

3
The first approach is empirical and the fatality rate is estimated using directly the ground motion intensity level and population exposed (Jaiswal et al., 2009). The casualty estimates using this approach are, in general, not satisfactory (Ranjbar et al. 2017), as shaking intensity is not directly linked to the number of deaths that also depends on vulnerability of the building stock. The second approach is semi-empirical and relates building damage to the number of casualties. In this work, casualty models based on the second approach will be employed. The casualty model proposed by So and Spence (2013) is a semi-empirical model and was tested against actual casualty data from previous events. The model is capable of considering different building classes (structural systems) to alter lethality rates for extensively damaged and collapsed buildings. In this casualty model, the number of people killed due collapse or extensive damage of buildings K for a given building class is defined as: where O is occupancy rate at the time of earthquake and P is the average number of people that normally reside in a building. d5 and d4 are the number of the buildings that collapse or are extensive damaged, respectively. L5 and L4 are the lethality rates representing the proportion of occupants killed for buildings that collapse ( d5 ) and extensively damaged ( d4 ), respectively.
The second model used in the analysis was proposed by Coburn et al. (1992), which considers only the collapse damage state for casualty estimations. Nonetheless, this model has additional parameters, such as occupants trapped by collapse and mortality post-collapse. The number of fatalities in this model K is expressed as follows: where D5 is the number of collapsed buildings, M1 represents the population per building, M2 is the occupancy rate at the time of earthquake, M3 is the number of occupants trapped by collapse, M4 and M5 are lethality rates for collapse and post-collapse, respectively.
The previous sections provided a general overview of the main components used in the proposed seismic risk assessment framework. To verify effectiveness and integrity of the framework, the following section uses the city of Adapazari in Turkey as a case study area.

Seismic hazard model for Adapazari
To assess the capability of the developed seismic risk assessment framework, the city of Adapazari in Marmara region (Turkey) was selected as a case study area. The Marmara region is located in one of the most seismically active regions in the world and was subjected to multiple earthquakes during the 20th century (Fig. 2). The North Anatolian Fault (NAF) lies across northern Turkey for more than 1500 km starting from Karliova in the east and extending to the Gulf of Saros in the west. Three main branches of NAF can be identified: the northern NAF (NNAF), central NAF (CNAF) and southern NAF (SNAF) branches. The NNAF dominates the tectonic regime of the Marmara Sea area.
The 1999 Kocaeli and Duzce earthquakes were the latest major events that ruptured the North Anatolian fault in the Marmara region. Due to its close proximity to the ruptured fault segment, Adapazari was one of the most severely damaged cities during the 1999 Kocaeli earthquake, suffering enormous economic losses, extensive structural damage and a high fatality rate. There is a high probability of another devastating seismic event to occur in the Marmara region in the foreseeable future , Murru et al. 2016. Despite this risk, there are no probabilistic seismic risk studies performed for the city. This highlights the need of the seismic risk assessment for Adapazari, as it is a vital resource for earthquake preparedness and risk mitigation (Erdik 2017).

Earthquake source zones
The earthquake source zone model is adopted from Sianko et al. (2020) and consists of background source zones (BSZs) and fault source zones (FSZs) that are used to generate synthetic earthquake catalogues. In total there are 17 BSZs (Fig. 3) and 25 FSZs ( Fig. 4) with unique earthquake recurrence parameters. BSZs are based on Gutenberg-Richter relationship, while FSZs are utilising characteristic magnitude M char . Aleatory  (Sianko et al. 2020) uncertainty is taken into account by randomizing key earthquake parameters with distribution functions during the generation of the synthetic catalogues.

Ground motion prediction equations (GMPEs)
Ground motion prediction equations have a big impact on both seismic hazard and seismic risk predictions . In this study, the GMPEs proposed by Akkar et al. (2014) and Boore et al. (2014) are used to predict ground motion intensities in terms of PGV, S a (T) , and S d (T) . The GMPE provided by Akkar et al. (2014) was based on data consisting of earthquake records mainly from Italy, Turkey and Greece, with the majority of records for strike-slip mechanism events with M > 7 obtained from earthquakes that occurred in the Marmara region. The GMPE provided by Boore et al. (2014) has a correction coefficient for different countries including Turkey. In Boore et al. (2014), inter-event and intra-event variabilities are determined for various magnitudes considering period, distance and soil conditions. Both GMPEs have a model for PGV, which is commonly used as ground motion intensity parameter in fragility functions for structures. Also both GMPEs are capable of considering various fault mechanisms and can utilise finite-fault ( R JB ) distance metric, which can be considered more accurate for earthquakes occurring on the faults than point-source ( R epi ). Local site effects can be estimated with relative ease as these GMPEs utilise widely available 30 m shear-wave velocities, V s30 . The V s30 is a widely used parameter to characterize seismic site conditions due to its relative availability and generally acceptable performance. Moreover, many recent GMPEs have site amplification functions based on V s30 values to take into account site conditions. The estimated topographic slope-based V s30 values are available globally through web-based open-access USGS map server. V s30 values obtained from topographic slope-based data can be used for large scale studies as shown by Riga et al. (2021). This is particularly important as a considerable part of the city of Adapazari is located in the areas prone to ground motion amplification. Epistemic uncertainty is addressed through employment of the logic tree, where 70% weight is given to and 30% weight to Boore et al. (2014). Higher weight is given to the GMPE developed by Akkar  (Sianko et al. 2020) et al. (2014) because the model contains a larger proportion of recordings from Turkey. Both GMPEs are capable of estimating S a (T) for a wide range of periods and satisfy criteria specified by Bommer et al. (2010) for use in PSHA.

Site conditions
Most of the area of the city of Adapazari is located on deep alluvial sediments deposited by the Sakarya River (Bray et al. 2004). The sub-surface soil is heterogeneous with big variations in soil layers. The soil generally consists of silty clays, silty sands, clean fine sands and gravels. The groundwater level varies seasonally, but on average it is 1-2 m below the ground surface. The shallow groundwater level contributed to the extensive occurrence of liquefaction in the city during the 1999 Kocaeli earthquake.
It is important for Adapazari to consider the site effects in the seismic hazard calculations, due to the presence of soft soil and possible site amplification effects. A firstorder approximation of the shear-wave velocities ( V s30 ) based on topographic slope model is employed in this work. V s30 map obtained from USGS shown in Fig. 5 is used in the PSHA performed for the city of Adapazari.

Exposure model for Adapazari
In this work, different sources of data are utilized to build the exposure model. The last published Building Census for Turkey was conducted in 2000. However, considering the rapid development of Adapazari since 2000, additional data is required. To address Fig. 5 Shear-wave velocities ( V s30 ) map for the city of Adapazari used in the PSHA this problem, aerial and satellite images are used to perform building footprint mapping to quantify the number and area of buildings in Adapazari. In addition, satellite images from past decades are examined to track the expansion of the city. This allows to approximately date the construction period of the structures in different parts of the city. The 2011 Population and Housing Census is used to determine average household size, average number of floors as well as proportion of buildings constructed before and after 1980.
In this work, a practical method based on building footprint mapping is used to find the number of buildings and their locations in Adapazari. The study area is divided in a grid consisting of 58 cells as shown in Fig. 6. Cell dimension is 0.01 × 0.01 degrees and all the data collected is at the cell level. The total number of buildings considered in the analysis is 47,283, where 31,067 were found from manual mapping (red coloured hatched area in Fig. 6) and 16,216 from assumed areas (green polygons in Fig. 6). In the assumed areas, the number of buildings is found by using density of buildings from adjacent mapped buildings areas. In this work, only residential buildings are used, thus public infrastructures, solely commercial or industrial buildings are not mapped or quantified.
Adapazari is the central district of Sakarya province for which the 2011 Population and Housing Census provides data about the building stock. From that census it is found that in the Sakarya province around 80% of the buildings are low rise (1-3 stories) and around 20% mid-rise (4-6 stories), with 2.5 being an average number of floors. 23% of the buildings were constructed before 1980, 60.7% were constructed after 1980 and 16.3% are of unknown date of construction. According to the 2000 Building Census, in Sakarya around 63% of the buildings were RC frames and around 35% bearing wall construction. The data from the Building Census are for a whole province including smaller town centers, where usually the higher proportion of buildings is bearing wall construction. Hence, it is safe to assume that the proportion of RC buildings in the city of Adapazari should be larger Fig. 6 Cells structure, mapped buildings and assumed polygons for the Adapazari in comparison to the province data. In this research it is assumed that 80% of the building stock is RC buildings and 20% are masonry with this ratio randomized within ± 10% using uniform distribution in the MC simulations.
The construction year of the buildings can be a useful parameter for selection of fragility curves, which in turn have an impact on seismic risk estimates. Seismic design Fig. 7 Satellite images of Adapazari and suburb areas in 1984 (a) and 2020 (b) guidelines in Turkey became more comprehensive after 1975 (Erdik et al. 2003), therefore in this study to reflect seismic resistance of the buildings, they are classified as pre-1980 (Low-Code) and post-1980 (High-Code). Satellite and aerial images from 1984, 1987, 1992, 1997, 2005 and 2020 are used to track development of the city and to assess the age of the buildings in each cell. Figure 7 shows an example of images with substantial time differences that were used for visual inspection. It can be noticed that the city expanded almost in all directions over the time period of 1984-2020. Figure 8 shows the result of visual inspection of satellite images represented with assigned code values, where 1 stands for high-code buildings and 2 is used for low-code buildings in the cell. This data will be used later in fragility functions capable of distinguishing low-code and high-code buildings.
Ground motion intensities from PSHA and scenario earthquake are calculated at the center of each cell. For simplicity, building locations inside the cells are represented by the centroids of the cells; a common assumption in seismic risk studies (e.g. Bommer et al. 2002, Erdik et al. 2003. The total number of buildings in each cell are given in Fig. 9. Remote street surveys are utilized at nodes of each cell to determine the average floor height and verify the age of construction estimated from satellite images. In total 40 buildings are surveyed per cell (10 buildings per node).

Economic value and population
Once the spatial distribution of buildings is determined, to perform earthquake risk assessment, the economic and population values need to be estimated. To find the roof area of buildings in the assumed areas (green polygons in Fig. 6) the distribution of the Fig. 8 Grid developed for the Adapazari with assigned code values, where 1 and 2 standing for high-code and low-code buildings respectively in the cell roof area across manually mapped buildings is utilised (Fig. 10). The total roof area is the sum of the roof area of the mapped buildings and the roof area estimated from the assumed areas for each cell. From the total roof area, total floor area is predicted by multiplying roof area by the average storey number in the cell. Finally, the cost of the buildings is found by multiplying the total floor area by average cost of construction per sq. meter. According to Turkish Revenue Administration (TRA 2021), the average cost of construction of RC residential buildings ranges from 1330 TL for first class buildings to 2130 TL for premium class buildings. The average cost of construction per sq. meter is 1083 TL for first class masonry buildings and is 1671 TL for premium class masonry buildings. In this study, the average cost of construction per sq. meter is assumed to be  1730 TL (230 USD) and 1377 TL (180 USD) for RC and masonry buildings, respectively (conversion rate 7.5 TL to 1 USD in March 2021). From the above calculations the total value of the building stock in the study area is estimated to be 5850 M USD.
According to the 2011 Population and Housing Census data, the average household size in Sakarya province was 3.9. Population per cell is calculated from the floor area in the cell multiplied by 70% (to exclude commercial buildings), then multiplied by the average household size and divided by the average dwelling area (130 sq. meters in Adapazari). The total population for the study area is found to be 538,178. This value is relatively close to the one provided on the Turkish statistics website (https:// www. tuik. gov. tr/). In 2020, the city had a population of 517,000 when Adapazarı city center (279,000), Serdivan (148,802) and Erenler (90,855) districts are combined to represent the study area in this paper.

Vulnerability model
The fragility and vulnerability models are the main sources of uncertainty in a seismic risk assessment procedure (Riga et al. 2017). Therefore, rigorous selection of these models is very important for accurate predictions of the seismic risk. In this study, seismic fragility models developed for Turkish building stock are considered.

Fragility curves
Most commonly used analytical fragility curves are not benchmarked against past earthquake events (Villar-Vega and Silva 2017), hence need further consideration before being used for this region. However, the fragility curves developed by Erdik et al. (2003) were obtained from post-earthquake damage data collected from Turkish building stock. These fragility curves are based on MSK-81 intensity and spectral displacement ( S d (T) ) for damage states slight, moderate, extensive and complete. The spectral displacement demand is calculated from spectral accelerations, using the displacement coefficient method given in FEMA 356 (2000). Furthermore, Erdik at al. (2003) proposed separate fragility curves for buildings constructed before and after 1980 (named pre-1980 and post-1980), which is useful for vulnerability model refinement in the areas where construction date is known. Moreover, these fragility curves include buildings construction type (RC and masonry) and buildings with various storey heights (low-rise, mid-rise, high-rise). Figure 11 shows an example of these fragility curves based on spectral displacement. Analytical fragility curves developed by Erberik (2008) are also used in this work. These curves are based on PGV and derived from 28 RC buildings extracted from a building database of approximately 500 buildings in Duzce, Turkey, which were constructed between 1962 and 1999. The fragility curves were verified against actual damage data. The examples of the fragility curves developed by Erberik (2008) are shown in Fig. 12.
Similar to Erberik (2008), Akkar et al. (2005) used a building database consisting of 32 low-rise and mid-rise typical RC buildings from Duzce (Turkey) to develop fragility curves based on PGV. Figure 13 shows the fragility curve for 3-storey RC buildings with slight, moderate and severe damage states. A comparison of Figs. 12 and 13 shows that the fragility curves developed by Erberik (2008) and Akkar et al. (2005) are quite similar and the main difference is observed at low PGV values for curves representing LS1, LS2 and light, moderate damage states, respectively. For example, the fragility curves developed by Akkar et al. (2005) predict no moderate damage for PGV values less than or equal to 30 cm/s, while fragility curves developed by Erberik (2008) predict a probability of exceedance of around 7% for LS2. It can be concluded that, in general, the fragility curves developed by Erberik (2008) predict higher probability of exceedance for moderate damage in comparison to those developed by Akkar et al. (2005). Kirçil and Polat (2006) developed analytical fragility curves based on peak ground acceleration (PGA), peak ground velocity (PGV) and pseudo-spectral acceleration (PSA) for RC frame buildings designed according to the Turkish Seismic Design Code published in 1975. They also proposed fragility curves for two different steel grades (S220 and S420). However, often there are no data available on the proportion of buildings constructed with each type of reinforcement. Hence, it may be more practical to

Fig. 12
Fragility curves for LS1-serviceability, LS2-damage control and LS3-collapse prevention limit states based on PGV for low-rise RC buildings (Erberik (2008) combine the fragility curves of two steel grades. One of the shortcomings of the fragility curves proposed by Kirçil and Polat (2006) is that only two limit states are considered: yielding and collapse. Moderate and extensive damage states are missing, and this makes the use of these fragility curves in seismic risk calculations difficult. Furthermore, low-rise buildings are only represented by fragility curves for 3-storey buildings, which may overestimate damage of 1-2 storey buildings. Figure 14 shows the combined material fragility curves for 3-storey RC buildings for yielding and collapse damage states.
In general, fragility curves based on PGA are not very accurate at short distances from the epicentre. At such distances, PGA values observed during the 1999 Kocaeli earthquake were under-predicted by most of the GMPEs due to the smoothness of  Kirçil and Polat (2006) for 3-storey RC buildings rupture and fairly low stress drop (Erdik 2001). On the other hand, the PGV estimates of GMPEs were similar to values observed in similar past events. This is particularly important as the case study area is located in close proximity to the fault that ruptured in that earthquake.

Consequence model and vulnerability curves
There are direct and indirect losses that can be caused by earthquake events. While direct losses are mostly associated with earthquake damage to the buildings, indirect losses can be caused by business and industry downtime due to recovery. In this study, only direct losses associated with structural damage are considered. Consequence models developed for Turkish building stock (e.g. Smyth et al. 2004, Crowley et al. 2005, Bal et al. 2008) are shown in Table 1 in addition to suggested values in HAZUS-MH MR1 (2003) for US building stock. It can be noticed that the consequence models for the Turkish building stock provide higher damage ratios than those proposed in HAZUS-MH MR1 (2003). According to the legal requirements in Turkey, buildings with extensive damages need to be demolished after the earthquake (Bal et al. 2008). Consequence models for Turkish building stock take this into consideration, by providing 100-105% of replacement cost for extensive damage,  while HAZUS considers 50% replacement cost for the same damage state. In the present study, the consequence model developed by Bal et al. (2008) is utilized as it also includes demolition cost for extensive and complete damage states. Figure 15 shows a comparison of vulnerability curves obtained with the convolution of fragility curves developed by Erdik et al. (2003) and consequence models proposed by Bal et al. (2008) and HAZUS-MH MR1 (2003). The effect of consequence models on the developed vulnerability curves can be observed from this figure. In Fig. 16 vulnerability curves for RC buildings with various number of storeys are obtained using Bal et al. (2008) consequence model and Akkar et al. (2005) fragility curves. The plateau for PGV values between 10 and 40 cm/s is due to the low contribution of moderate and severe damage fragility curves to MDR for this PGV range.

Seismic risk analysis results for Adapazari
MC-based probabilistic seismic risk analysis is performed for the city of Adapazari and results are presented in this section. In addition, a scenario earthquake similar to 1999 Kocaeli earthquake is modelled and a seismic risk analysis is performed using a scenario earthquake to compare results with the observed damage and casualties.

Seismic hazard
An event-based PSHA is performed for the case study area using background and faults source zones as explained in Sect. 3.1. Ground motion intensities (such as PGA, PGV and S_a(T)) are calculated from GMPEs considering soil amplification with intra-event and inter-event variabilities. Hazard maps for Adapazari are developed for probability of exceedances 50%, 10% and 2% in 50 years (return periods of 72, 475 and 2475 years, Fig. 16 Comparison of vulnerability curves obtained by convolving the consequence model developed by Bal et al. (2008) with the fragility curves developed by Akkar et al. (2005) for ordinary 3,4,5-storey RC buildings  Figures 17 and 18 present PGA and PGV values obtained from hazard analysis at the centre of each cell for the specified return periods. The results of this work show that the PGA values for Adapazari range between 0.65 and 0.83 g with an average value of 0.74 g, for a 475-year return period; and between 1.06 and 1.44 g with an average value of 1.23 g, for a 2475-year return period. In the seismic hazard study performed by Erdik et al. (2004) for the Marmara region, PGA values for Adapazari were predicted to be in a range between 0.6 and 0.8 g and in a range between 1.0 and 1.5 g for 475 and 2475 years return periods, respectively. These values are in line with PGA values predicted in this work. In a more recent study performed for the eastern part of Marmara region, Gülerce and Ocak (2013) predicted PGA values for Adapazari to be between 0.6-0.8 g and 1.0-1.2 g for 475 and 2475-year return periods, respectively. The PGV values for Adapazari range between 61 and 89 cm/s with an average value of 79 cm/s, for a 475-year return period; and the PGV values range between 100 and 141 cm/s with an average value of 128 cm/s, for 2475year return period.

Seismic risk
Event-based probabilistic seismic risk analysis results are presented in this section. Figure 19 shows the loss curves for the city of Adapazari in terms of MDR and economic loss in USD using Erdik et al. (2003) and Erberik (2008) fragility curves. Building period was randomised when Erdik et al. (2003) fragility curves were used. Uncertainty of MDR is not considered as it slightly affects seismic risk analysis results. It can be seen from the loss curves that both models result in similar MDR values for smaller return periods. The loss curves based on fragility curves developed by Erberik (2008) predict higher economic loss for larger return periods. In the development of the latter curves, all buildings in the study area were assumed to be RC. Figures 20 and 21 show MDR maps using the fragility curves developed by Erdik et al. (2003) and Erberik (2008) for return periods of 72, 475 and 2475 years, respectively. As expected, the MDRs based on fragility curves developed by Erdik et al. (2003) in comparison to those developed by Erberik (2008) result in lower and higher values in a recently constructed and old areas of the city, respectively. Figure 22 presents a lethality map for Adapazari based on fragility curves developed by Erdik et al. (2003) and casualty model proposed by So and Spence (2013) for various Fig. 19 Loss curves for Adapazari based on fragility curves developed by Erdik et al. (2003) and Erberik (2008) 1 3 return periods. Occupancy rate in the casualty model is assumed to be 0.9 (So and Spence 2013), which is typical value for night time. Table 2 shows total lethality predictions for the city of Adapazarı based on casualty models proposed by So and Spence (2013) and Coburn et al. (1992) for return periods of 72, 475, 2475 years. These results show that lethality rates for night time earthquakes are 3 to 4 times higher than those for day time earthquakes. Coburn et al. (1992) model results in higher lethality rates for all return periods in comparison to So and Spence (2013) model.

Scenario earthquake and comparison with Kocaeli earthquake
A scenario earthquake with characteristics similar to 1999 Kocaeli earthquake is generated using the hazard analysis tool developed in a previous study (Sianko et al. 2020). The fault rupture model is represented by a single rectangular plane with magnitude of M w = 7.4 . The GMPEs mentioned in Sect. 3.2 are used for the generation of the scenario earthquake.  Table 2 Lethality estimates based on So and Spence (2013) and Coburn et al. (1992)  To verify the developed framework, seismic risk assessment is performed for Adapazari using the generated scenario earthquake. The importance of scenario earthquake risk assessment is to predict the damage distribution of the buildings' portfolio caused by the modelled scenario earthquake (Kohrangi et al. 2021b). Figure 23 shows the spatial distribution of median PGA values obtained from the scenario earthquake. The PGA ranges from 0.31 g in the northern part of the city and gradually increasing to 0.52 g towards the fault in the south. The PGA recorded during 1999 Kocaeli earthquake at Sakarya station was 0.41 g, but according to Kudo et al. (2002) the station is located on very stiff soil and ground intensities in the alluvial basin, at the central part of Adapazari, could be substantially different from those recorded by the station. As a result, the station does not represent the site conditions at the city center, where the most of the damage was observed (Bakir et al. 2002). In the hazard analysis, the PGA values in the city center of Adapazari is estimated to be in the order of 0.3-0.4 g as suggested by previous studies (e.g. Sancio et al. 2002;Yakut et al. 2005). Following the 1999 Kocaeli earthquake, the municipality of the Adapazari performed a damage assessment for the buildings in the city . The total number of buildings assessed by the municipality was 23,914 and these buildings were categorised into two damage groups. The first group represents buildings that can be repaired (light and moderate damage) and the second group is for buildings that are required to be demolished (extensive and complete damage). Table 3 provides a comparison of percentages of buildings with extensive and complete damage at district level obtained from the scenario earthquake and survey performed by Adapazari municipality following the 1999 Kocaeli earthquake. It can be observed from this table that the results of the scenario earthquake hugely overestimate damage for the districts located in the south of the study area (e.g. districts 1-4). This is an unexpected finding as those districts are located closer to the fault rupture  Mollamahmutoglu et al. (2003) and GMPEs are predicting higher ground motion intensities in these areas considering the soil conditions. Therefore, observed damage in these areas is expected to be higher. On the other hand, the seismic risk analysis performed using the scenario earthquake underestimates the observed damage in the central districts (e.g. districts 20-22). This can be partially explained by the liquefaction phenomena observed in these areas (Fig. 23). A large proportion of buildings located in these areas suffered extensive and complete damage due to liquefaction instead of strong ground motion itself. In addition, according to Mollamahmutoglu et al. (2003) even buildings designed according to the seismic regulations suffered from severe displacements. Hence, poor predictions for central districts might be due to GMPEs not being able to predict soil amplification correctly as occurred in 1999 Kocaeli earthquake. Nevertheless, Bakır et al. (2005) predicted higher ground motion intensity for the north and central parts of the city using microzonation and soil response analysis than presented in the current study with GMPEs. Disagreement in results can be also due to seismic risk procedures do not consider damage due to liquefaction. It can be concluded that it is not possible to predict damage patterns observed during the 1999 Kocaeli earthquake using conventional methods that are developed for large scale studies. Similar conclusion for scenario earthquakes for Adapazari was drawn by Yakut et al. (2005).  Table 4 shows overall comparison between scenario earthquake and survey data in terms of % of total number of buildings. As can be seen, predictions for extensive and complete damage states obtained from the proposed risk analysis for the scenario earthquake correlate well with the damage data collected following the 1999 Kocaeli earthquake. On the other hand, the number of buildings with slight and moderate damage states are highly overestimated by the proposed risk analysis method in comparison to the survey data. During the 1999 Kocaeli earthquake, numerous surface manifestations of liquefaction, in the form of sand boils and lateral spreading, were observed in Adapazari. The softened and/or liquefied soil acted as an isolator dissipating the energy at the foundation level and reduced shaking damage to the buildings (Bakir et al. 2002). Therefore, most of the buildings with slight and moderate damages predicted by the earthquake risk assessment procedure did not suffer any damage during the 1999 Kocaeli earthquake. This highlights an urgent need for the integration of the effect of liquefaction susceptible soil conditions for buildings located on such soils into the risk analysis. Figure 24 shows the damage distribution in Adapazari obtained from risk analysis using the fragility curves developed by Erdik et al. (2003). Complete and extensive damage states are mainly observed in the areas where the density of the building population is high.
The lethality predictions for the scenario earthquake estimates 3817 and 5298 deaths using So and Spence (2013) and Coburn et al. (1992) casualty models, respectively. The occupancy rate at the time of the earthquake is assumed to be 0.9 (as the Kocaeli earthquake occurred at night time). Bar-Dayan et al. (2000) and Margalit et al. (2002) stated that the 1999 Kocaeli earthquake killed 2680 people and approximately 5300 injuries in Adapazari. However, Bakir et al. (2002) reported fatalities as 3684 for the same earthquake. The population of Adapazari was 283,752 according to the census taken in the 2000s. In this work, the population of Adapazari is estimated as 517,000. It can be concluded that both casualty models used in this work give reasonable results.

Conclusions
In this work, a probabilistic seismic risk assessment framework based on MC simulations is developed. The developed framework is capable of producing seismic risk maps and loss curves. The uncertainty of input parameters in seismic risk is treated by employing logic tree and randomization processes in MC simulations. An exposure model for the study area is obtained with building footprints mapped from aerial and satellite images and remote street survey. The developed earthquake risk assessment framework is applied to the city of Adapazari in Turkey using mostly easily accessible data. The results for Adapazari are presented with PGA and PGV hazard maps, MDR distributions for various return periods and loss curves. The developed risk maps for the city of Adapazari can help stakeholders and decision-makers to reduce future earthquake related losses.
One of the biggest challenges in the development of a reliable seismic risk framework is the verification of obtained results against post-earthquake data. The developed framework was verified by comparing the building damage predicted for a scenario earthquake with those observed after the 1999 Kocaeli earthquake. In addition, casualty models were employed in the earthquake risk analysis for the estimation of fatalities in the city of Adapazari due to the scenario earthquake. Risk analysis results demonstrated that the developed earthquake risk assessment framework can predict extensive and complete damage states in Adapazari very reasonably for the scenario event in comparison to damage observed after the real event. However, slight and moderate damage states are overestimated by the developed procedure. This could be partially explained with the soil conditions of Adapazari, which lies over soft and liquefiable silts and sands. During the 1999 Kocaeli earthquake, the liquefied soil could have dissipated some energy acting as an isolator and reduced the shaking damage to the buildings.
To conclude, the developed framework can serve as an efficient tool for the assessment of seismic risk, but it should be used with caution in areas with complex geological conditions that are prone to liquefaction. There is an urgent need for the development of a Probabilistic Liquefaction Hazard Analysis (PLHA) procedure based on MC simulations to derive liquefaction hazard curves and PLHA maps for seismic prone regions. Using the PLHA procedure, soft soil conditions and liquefaction effects can be integrated into the probabilistic seismic risk analysis for buildings located on such soil conditions. Also, a first-order approximation of the shear-wave velocities ( V s30 ) based on topographic slope model is employed in this work. In a future work, a more accurate estimation of basin effect will be included in the PSHA procedure.
Author contributions All authors contributed to the study conception and design. Data collection was performed by IS and ZO. The development of codes and analysis were performed by IS. Outputs were checked by ZO, IH, and KP. The first draft of the manuscript was written by IS and ZO and all authors commented on previous versions of the manuscript. All authors read and approved the final manuscript.