Significant Seismic Risk Potential From Buried Faults Beneath Almaty City, Kazakhstan, Revealed From High‐Resolution Satellite DEMs

Major faults of the Tien Shan, Central Asia, have long repeat times, but fail in large ( Mw 7+) earthquakes. In addition, there may be smaller, buried faults off the major faults which are not properly characterized or even recognized as active. These all pose hazard to cities along the mountain range front such as Almaty, Kazakhstan. Here, we explore the seismic hazard and risk for Almaty from specific earthquake scenarios. We run three historical‐based earthquake scenarios (1887 Verny Mw 7.3, 1889 Chilik Mw 8.0 and 1911 Chon‐Kemin Mw 8.0) on the current population and four hypothetical scenarios for near‐field faulting. By making high‐resolution Digital Elevation Models (DEMs) from SPOT and Pleiades stereo optical satellite imagery, we identify fault splays near and under Almaty. We assess the feasibility of using DEMs to estimate city building heights, aiming to better constrain future exposure datasets. Both Pleiades and SPOT‐derived DEMs find accurate building heights of the majority of sampled buildings within error; Pleiades tri‐stereo estimates 80% of 15 building heights within one sigma and has the smallest average percentage difference to field‐measured heights (14%). A moderately sized Mw 6.5 earthquake rupture occurring on a blind thrust fault, under folding north of Almaty is the most damaging scenario explored here due to the modeled fault stretching under Almaty, with estimated 12,300 ± 5,000 completely damaged buildings, 4,100 ± 3,500 fatalities and an economic cost of 4,700 ± 2,700 Million US dollars (one sigma uncertainty). This highlights the importance of characterizing location, extent, geometry, and activity of small faults beneath cities.

• Digital elevation models derived from high-resolution satellite imagery can map active faulting near cities and determine building heights • Scenario risk calculations show a moderate earthquake on a fault in north Almaty would cause considerable damage and loss due to proximity • Properly characterizing fault location and geometry close to cities is key to quantifying the relative level of seismic hazard and risk

Supporting Information:
Supporting Information may be found in the online version of this article.  (Farr et al., 2007). Population estimates from United Nations, World Urbanization Prospects 2018 revision (United Nations, 2018); in 2020 the population is  E 1,896,000. (b) Tectonic setting of Almaty in the northern Tien Shan. GPS velocities relative to stable Eurasia with 95% confidence ellipses, from Zubovich et al. (2010). Active faults from Mohadjer et al. (2016) in red and Grützner et al. (2017) in light orange; locations of historical earthquakes in white from (Kalmetyeva et al., 2009); yellow line shows the boundary between Kyrgyzstan and Kazakhstan. (c) RGB Landsat image (2017) of the city of Almaty, active faults in light orange from Grützner et al. (2017) and fault to the north of Almaty mapped in this study shown in dark orange.
AMEY ET AL.

10.1029/2021EA001664
3 of 25 within living memory. Additionally, many large cities in this region formed near or on the mountain range fronts, between the plains and high mountains, where the mountains offer protection. These settlements have grown in the past century into major cities and have expanded toward and onto the mountain range fronts and consequently are now very exposed and vulnerable to earthquakes (Figure 1). This urbanization additionally results in modification of the landscape and potentially obscures or obliterates the geomorphic markers of faulting activity, thus making it harder to characterize the hazard.
In this study, we explore damage and losses due to earthquake scenarios for Almaty, the largest city in Kazakhstan, formerly known as Verny (Figure 1). It offers a particularly poignant study site as it has been leveled several times by large earthquakes, including the 1887 Verny earthquake in which all stone buildings were reported to be destroyed (Hay, 1888), with 180 fatalities in Verny town and 328 fatalities in total. At the time, the population was 25,000, whereas due to rapid growth and urbanization the population is now ∼1.9 million. The city is built on an alluvial fan and expansion extends further northwards into the Kazakh plain, southwards onto the range front and also outwards, subsuming other towns along the east-west striking range front (Figure 1a). The evidence for active faulting has been laid out in Grützner et al. (2017), who found fault scarps of 2 to  E 20 m height along the Zailisky rangefront East of Almaty, suggesting the faults have been active and ruptured to the surface since the Last Glacial Maximum.
Here we create high resolution Digital Elevation Models (DEMs) using optical stereo satellite imagery (Pleiades and SPOT) to identify topographic scarps of recent ruptures that may be present in the city, but are now hidden due to the urban expansion. We aim to better understand which potentially active faults may pose a risk to the city. Additionally, we show a proof of concept that building heights in Almaty can be remotely calculated from high resolution DEMs, which in the future may improve exposure models across a city where such information may be hard to gather. We use a deterministic seismic hazard modeling approach to assess the degree of ground shaking hazard and likely range of damage and loss, due to plausible earthquake scenarios. We use the Global Earthquake Model's (GEM) OpenQuake engine (GEM, 2018;Pagani et al., 2014;Silva et al., 2020) and run seven earthquake scenarios; four are hypothetical scenarios based on our fault mapping and three are historical earthquakes known to have caused significant damage to Almaty. The historical earthquake scenarios provide a baseline to understand the impact if these past earthquakes were to happen today. We extend the work of Mosca et al. (2019) by calculating not only seismic hazard in terms of estimated levels of ground shaking, but additionally the potential risk to Almaty in terms of damaged buildings and loss of life. We present the damage to residential buildings and the human and economic losses due to these earthquake scenarios, and identify which scenario is most damaging.

Identifying Fault Scarps and Building Heights From High Resolution DEMs and Field GPS Measurements
Freely available DEMs enable us to map the major active faults along the range front at a moderate regional scale based upon their geomorphic expression in the landscape (Grützner et al., 2017), including SRTM at 30 and 90 m resolution, the ASTER Global 30 m DEM and ALOS world DEM (AW3D, 30 m). These are not of a high enough resolution to be able to identify individual buildings or subtle geomorphic expressions of fault splays within the city that could potentially cause damaging earthquakes ( Figure 2). Instead, we use the higher resolution 12 m TanDEM-X DEM (Borla Tridon et al., 2013) and create higher resolution DEMs from SPOT and Pleiades stereo imagery.
To derive topography, we process along-track stereo SPOT, two stereo Pleiades scenes and one tri-stereo Pleiades set of panchromatic optical images (all supplied with Rational Polynomial Coefficients, RPCs). We use Agisoft Metashape version 1.6.2 to build dense clouds and create DEMs for the Pleiades images, after initially using ERDAS IMAGINE Photogrammetry software (ERDAS IMAGINE 2016), but found tall buildings were not resolved using ERDAS. As part of this workflow, tie-points are automatically identified between images, from which triangulation is calculated. From this, a point-cloud of points with latitude, longitude, and elevation is produced. We process this to Ultra High resolution (which uses the native resolution of the input images) and use "mild" depth filtering, as we found that increased filtering removes tall buildings as anomalous noise. We used ERDAS software standard photogrammetry workflow to process the AMEY ET AL.
10.1029/2021EA001664 4 of 25 SPOT images; we found ERDAS produced better results than Agisoft for SPOT, discussed in Section 4. We note that AMES stereo pipeline is a freely available software not used here, as it was not designed for urban environments. Due to potential geolocation errors between the satellite-derived DEMs, we coregistered each DEM to the 12 m Tandem-X DEM following the geolocation procedure of (Nuth & Kääb, 2011). The coregistration shifts are reported in Table S1 and all DEMs are referenced to the ellipsoid. Pointclouds of the DEMs derived in this study are available to download on OpenTopography (Amey, 2021a;Amey, 2021b;Amey, 2021c).
We remove buildings from the pointcloud to investigate whether there are any subtle splays within the city, by filtering the buildings using LAStools (Isenburg, 2020), script availability listed in acknowledgments. This is done by first removing high and low outliers in 3.5 m cells, then classifying the lowest non-noise points, and classify the ground in steps of 25 m (as is recommended for a city [Isenburg, 2020]). We output the digital surface model (DSM, including buildings) and digital terrain model (DTM, ground level) at 1.5 m resolution for Pleiades and 4 m for SPOT. Digital Elevation Models (DEMs) over the same part of Almaty city, where each pair of panels has colored topography and a zoom-in of a  E 2 km region of the eastern city in gray hillshade. Left column are freely available DEMs, including: Shuttle Radar Topography Mission (SRTM) at 90 m resolution (Reuter et al., 2007), SRTM at 30 m resolution (Farr et al., 2007), ASTER Global DEM at 30 m (Tachikawa et al., 2011), the AW3D (ALOS World 3D) at 30 m (Takaku et al., 2016, original data are provided by JAXA). Right column are commercial DEMs including TanDEM-X at 12 m (Borla Tridon et al., 2013), and DEMs that we have derived from SPOT and Pleiades stereo optical imagery at 4 and 1.5 m respectively. The lower resolution DEMs are able to pick out large, mountain-range scale fault structures but are unable to measure houses, whereas the higher resolution clearly picks out the outline of buildings. Note that the Pleiades tri-stereo scene shown here does not have coverage as far south as the other DEMs.
10.1029/2021EA001664 5 of 25 Grützner et al. (2017) previously identified a fault scarp within Almaty from field observations and analysis of SRTM and SPOT DEMs. The topographic step is consistent on parallel roads up to 1 km apart, suggesting this is the trace of an active fault striking east-west, consistent with a fault stepping into the basin from the main mountain range ( Figure 6). Here we investigate this further with Pleiades stereo and tri-stereo DEMs and ground-truth with Global Navigation Satellite System field measurements (GNSS), in order to determine along-strike extent and magnitude of offset.
We use a Leica GS18 GNSS attached to a Leica CS20 hand-unit, with the Smartlink function in order to use the GS18 as a rover without having a base-station for real-time kinematic operation. Corrections are supplied via geostationary augmentation satellites, to which the link must be maintained or positional accuracy will drop dramatically. As such, if connection was lost we maintained a stationary position till the link was re-established. The GS18 was mounted on a 2 m pole and took a measurement every second. We position the GS18 through the window with the antenna above the car, and measured the height off the ground to correct for this offset. We drove along north-south roads shown in Figure 3 multiple times, in which there were noticeable offsets. In Figure 3 we only plot the pass with the lowest uncertainties on the datapoints, that is, for passes in which the satellite signal was not lost, which otherwise leads to extremely high (up to 9 m) uncertainties and wildly varying positions. This means in profile P1 we find points greater than 3.5 km distance along the profile, as this is where the road went through an underpass, which caused the GNSS to lose satellite connection leading to a lack of reliable data. We also removed points with uncertainties greater than 4 m in profile P2 (at more than 3.8 km along the profile), and those with uncertainties greater than 2.5 m in profile P3 (at less than  E 1.1 km along profile). The mean error of the four GNSS profiles are 1.2, 1.7, 2.2, and 1.5 m respectively. We also took profiles along these roads in the satellite-derived DTMs, taken with bins of 5 m and swaths of 2 m for Pleiades, and bins of 10 m and swaths of 4 m for SPOT.
Profiles P1 to P4 show subtle but measurable offsets which we attribute to the surface trace of a fault splay under Almaty (Figure 3), the magnitude of which we have calculated separately for each data set (Table 1). We note in some profiles there is an offset between the GNSS and DTM elevations-both datasets are referenced to the ellipsoid and we have not done any post-processing/co-registration of the GNSS here. We also estimate the offset of the fault for six profiles for which we do not have GNSS ground-truthing ( Figure S1). We attribute the differences in magnitude of the offset between satellite DTMs for the same profile, to how well the buildings have been removed from the DEMs in the lastools processing as well as errors in the satellite-derived DEM. Using the location of these offsets, we map out the fault originally identified by Grützner et al. (2017) through Almaty for hazard modeling.
In addition, (Grützner et al., 2017) identified folding in the north of Almaty city, partially in the city's Alatau district. Using a TanDEM-X DEM we map the extent of this fold (Figure 4). We suggest that this it is a fault-propagation fold, caused by a fault dipping to the south that extends under Almaty. The city building grid system runs approximately North-South and we take profiles along the major roads. The map is roughly the extent of the Pleiades tri-stereo and outlines show the extent of the two stereo (green and orange) Digital Elevation Models (DEMs). (b) to (e) Show profiles in the city from SPOT (blue), Pleiades stereo (light green and orange) and Pleiades tri-stereo (red) with field-measured Global Navigation Satellite System (GNSS) (black with error bars in gray). The best-fit lines to the GNSS data are shown in dark green, with the estimated fault offset from GNSS measurements (see Table 1.) Note that in panel (a) we plot the DEM, but the profiles have been taken on the digital terrain model, that is, the terrain model in which buildings have been removed from the DEM.

Exploring Earthquake Scenarios for Almaty With OpenQuake
OpenQuake Engine is an open-source software tool for earthquake hazard and risk modeling Silva et al., 2014). We use "scenario" mode to model specific, plausible earthquakes for Almaty based upon past and credible potential seismic events.
The OpenQuake engine can be considered in two components: hazard and risk. For scenario hazard calculations, an input earthquake rupture is used to calculate a distribution of ground shaking in a given region. The earthquake is described using a finite fault plane, with a given magnitude, rake and geometry. Ground motion is estimated using ground motion prediction equations (GMPEs): empirical formulae that predict intensity measures as a function of distance to the fault. Here we calculate peak ground acceleration (PGA), the maximum acceleration that occurs at a given location, and spectral acceleration (SA), the maximum acceleration that occurs at a given location on a 1-D damped harmonic oscillator for a given period. Peak ground acceleration gives the movement felt by a particle on the ground; spectral acceleration approximates acceleration experienced by a building. Here we use three GMPEs developed for shallow crustal earthquakes, see details in the original papers Campbell & Bozorgnia, 2014;Chiou & Youngs, 2014). The scaling of ground motion (PGA) with distance is shown in Figure S2. These GMPEs use the distance to surface rupture ( rup E R ) in their calculations, meaning distance from an exposure point is calculated using closest distance to the fault plane, whether at the surface (for the foot-wall) or at depth (for the hanging-wall). Exposure points in the hanging wall experience elevated shaking that decays with the distance from the surface to the fault at depth. These three GMPEs characterize the hanging-wall effects on results of 1-D finite-fault kinematic simulations from Donahue and Abrahamson (2014), though the exact scaling is calculated differently in each case. Choice of epicenter is irrelevant to these calculation as it is not used in these GMPEs. Rake is used to determine the style of faulting-either reverse, normal or strike-slip, and so in selecting the rake for hypothetical scenarios we choose purely reverse faulting, except for the deep earthquake in the Kazakh platform. For the earthquake scenario north of Almaty in the Kazakh platform, as this is a thick, stable crustal region we use the Somerville (2009) GMPE that is developed for stable cratonic regions (central and western Australia).
GMPEs use time-averaged shear-wave velocity in the top 30 m ("Vs30") to estimate the near-surface amplification of ground shaking caused by soil and unconsolidated sediment. We use the global topographic slope derived Vs30 map from Allen and Wald (2009), with a resolution of 30 arc-seconds (∼1 km), in the absence of measured Vs30 in Almaty. This assumes that steep slopes are formed of high-velocity, stiff materials and that basins are formed of low-velocity, soft sediments. For each exposure point, at which we calculate hazard and risk, the nearest Vs30 value is used. Figure S2 shows the distribution of exposure points (1-1.5 km in the center of Almaty, and up to 11 km apart in the south), with the closest Vs30 point to an exposure point giving a good approximation of the whole area. In theory with sufficient information available (e.g., shear wave velocity profiles and local geology) we could either compute directly a map of Vs30 or recalibrate the Vs30-slope relationship. The first option would be in any case the preferred one since inferring Vs30 from Pleiades tristereo (m) Pleiades stereo (m) Pleiades stereo (m) SPOT (m)  Figure 3), Using GNSS, SPOT, Two Pleiades Stereo Scenes and One Pleiades Tri-Stereo Scene 7 of 25 the slope gives only a rough approximation, and gains from deriving a higher resolution map would be limited. We acknowledge the limitations of using Vs30 derived from slope, especially for local risk assessments (Lemoine et al., 2012;Narciso et al., 2013), but in absence of more precise information we consider this the best option available. The Allen and Wald (2009) study also does not include analysis on whether this is a proxy for dry or water-saturated soil, and as such also does not account for any seasonal changes in site effects as has been suggested is important for Almaty (Alshembari et al., 2019).
GMPEs contain uncertainties and so to sample this, we use a Monte-Carlo sampling approach to randomly sample 1,000 realizations. We sample the logarithm of the ground motion using a Gaussian distribution with truncation level of three standard deviations. We assume that the ground motion should be correlated, and use the spatial correlation of Jayaram and Baker (2009)  From these intensity measures (PGA or SA) the OpenQuake engine risk module then calculates damage and losses. This requires an exposure model, which describes in detail the assets being exposed to the hazard; at given locations the building class, number of stories, number of buildings, cost of replacement and number of occupants for daytime and nighttime are listed. For each specific building type, OpenQuake calculates damage using corresponding fragility curves and losses (cost and fatalities) using corresponding vulnerability curves. We use fragility and vulnerability curves formulated for Central Asian building types (Martins & Silva, 2020). (2020) using screening surveys from local structural engineers, data from the population and housing censuses and high-resolution satellite images. The exposure is spaced every 1-1.5 km in the center of Almaty ( Figure S2). We use only the residential exposure data, with 2015 census data (details in Pittore et al., 2020). This exposure model potentially over-estimates the number of high-rise buildings, which we explore further in the discussion. The main types of building in Almaty are shown in Figure 5. For each district, the majority of buildings are unreinforced masonry, followed by reinforced masonry; very few of the total building stock are wooden or steel constructions. The Medeu district is largest and has the most residential buildings (Table S2). As the district of Medeu is large and has a distinct split in building density due to the mountain range, we split this district into two regions following the steep change in topography ( Figure 6).
Seismic hazard scenarios for three of the largest historical earthquakes have been calculated by Mosca et al. (2019). Here we modify their input ruptures for these earthquakes, suggest further earthquakes that may present hazard and additionally calculate the losses and damages in Almaty due to these earthquakes. Mosca et al. (2019) selected GMPEs that implement a mixture of the rup E R (Chiou & Youngs, 2014) and JB E R rupture-site metrics (Akkar et al., 2014;Boore et al., 2014).
R is the Joyner-Boore distance, which is the distance between the site and the closest point of surface projection of the rupture surface, that is, all sites located vertically above the rupture surface have JB E R equal to 0 and the largest values of shaking; the depth of the fault below the surface is not taken into account. Here, we use the GMPEs that use rup E R , Abrahamson et al. (2014); Chiou and Youngs (2014); Campbell and Bozorgnia (2014) due to the complex ruptures we model, as we find this more appropriate for faults with a gentle dip, such as the ones in proximity to Almaty.

Earthquake Scenario Selection
We run seven earthquake scenarios, described in detail below and shown in Figure 6, with full parameters given in Table S2. This suite of scenarios is comprised of, first, three historical earthquakes which we model to demonstrate the impact these earthquakes would have on present-day Almaty, given its current population and the distribution and typology of buildings. The other four we base on potential earthquakes we believe to be a plausible risk to Almaty. This is based on past work (Grützner et al., 2017), and the geomorphic assessment of active faults close to Almaty done in this study. Additionally we explore the effect of a large but deep earthquake further north, in the Kazakh platform (based on seismicity from Sloan et al., 2011).

M w 7.3 Verny 1887
Almaty (originally known as Verny) was founded in 1854, and suffered a major earthquake on June 9, 1887. The population at the time was  E 25,000 and the town consisted of 1,839 brick houses and 840 wooden ones, though less than 30 were two-storied (Hay, 1888). This earthquake killed 180 people in Verny and 328 in total, and all stone buildings were reported to be destroyed (Hay, 1888). No surface rupture was reported or has since been identified; it is not known whether this is because the fault slip did not reach the surface, or if it ruptured high in the mountain range and has been subsequently destroyed by erosion. Here we seek to quantify how much greater the damage and fatalities would be if this earthquake were to happen today, due to its magnitude and likely proximity.
We model the rupture as M w E 7.3 (Kondorskaya et al., 1982in Tatevossian, 2007 with a 82 km length fault using scaling relations from Wells and Coppersmith (1994), with the fault extending from the surface to 20 km depth, a dip of 45° south, rake of 90°, and strike of 65° based on isoseismals from Mushketov (1891) in Tatevossian (2007). This puts the fault at the southern end of Almaty city, intersecting the districts of Medeu, Bostandyk, and Nauryzbay ( Figure 6).

M w 8.0-8.3 Chilik 1889
Two years after the Verny earthquake, Almaty was damaged further by the July 11, 1889 Chilik earthquake. This is potentially one of the largest instrumentally recorded continental earthquakes, with an estimated magnitude of  E 8 (K. E. Abdrakhmatov et al., 2016;Bindi et al., 2014;Krüger et al., 2015Krüger et al., , 2018. Of the three historical earthquakes we present, this is suggested to be the largest magnitude, but caused the least damage  (Pittore et al., 2020). Because the Medeu district stretches into the mountains, later in this paper we have split the Medeu district into Medeu north and Medeu south, following the step in topography (shown by a white dashed line). (b) The predominant building types in Almaty in each district. In every district, the principal building material is unreinforced masonry, which typically performs poorly during earthquake shaking (high fragility). Each district has a high proportion of reinforced masonry structures, followed by reinforced concrete and finally confined masonry. There are very few wooden or steel structures (less than 1% of total).  from Mushketov, 1891. This may be due to its further distance from Almaty (∼100 km southeast) and related to a change in building type in Almaty, following the 1887 earthquake.
We use surface ruptures mapped by K. E. Abdrakhmatov et al. (2016) who extended the original field mapping by Crosby et al. (2007), Tibaldi (1998), and Tibaldi et al. (1997) The surface ruptures show a complicated rupture pattern, with two right-lateral strike-slip faults dipping at 90° (Beshkaragai and Kurmentey rupture) joined by a left-lateral strike-slip fault (Saty rupture) dipping at 40° N (K. E. Abdrakhmatov et al., 2016). We simplify the mapped rupture pattern and model the earthquake as purely right-lateral (rake = 180°), with a magnitude of M w E 8.0 (rather than 8.3, as is favored by Krüger et al. (2015) from amplitude comparison of long period surface waves) and rupture stretching from the surface to 20 km depth. The strikes of the three ruptures are shown in Figure 6; the two vertical faults have a strike of approximately 100° and 110°, the fault dipping at 40° N has a strike of approximately 240°.

M w 8.0 Chon-Kemin 1911
The Chon-Kemin earthquake on January 3, 1911 is reported to have killed 452 and injured 740 (Delvaux et al., 2001). Following original mapping of the ∼150 km long rupture by Bogdanovich et al. (1914), this has been extensively and systematically mapped in detail by Arrowsmith et al. (2017), and has been documented to have caused many landslides (Delvaux et al., 2001). This earthquake resulted in significant soil Figure 6. Earthquake scenarios used as inputs for OpenQuake in this study, where red is the surface rupture, blue shows the down-dip extent of the faults and the black arrow indicates slip vector. The Verny, Chilik, and Chon-Kemin scenarios are historical earthquakes known to have caused damage to Almaty. We also suggest four reasoned earthquake scenarios, based largely upon the existing knowledge of presence and location of active faults in and near the vicinity of the city. These include a deep, large earthquake in the Kazakh platform north of Almaty and three smaller earthquakes on identified fault splays close to Almaty (districts of which are denoted by the black outlines). All earthquake scenarios are suggested to extend from the surface to 20 km depth, except the Kazakh platform earthquake, which is buried, and extends from 30 to 50 km depth.
AMEY ET AL.

10.1029/2021EA001664
10 of 25 deformation and ground failure within Almaty itself, and is suggested to be because a shallow frozen soil layer induced liquefaction at depth (Alshembari et al., 2019).
As this fault rupture has been well mapped (Arrowsmith et al., 2017), we use this earthquake to examine the effect of our knowledge of the degree of fault complexity on the estimated range of losses. In the discussion, (Section 5.7) we run three different scenarios, parameterizing the rupture as one, two or a complex fault rupture.
To compare with the other scenarios, we use the two-fault parameterization for the Chon-Kemin earthquake as a complex fault rupture may not necessary reflect the complexity at depth. We define two fault planes that fit the majority of Arrowsmith et al. (2017)'s mapped surface ruptures. We use a magnitude of 8.0, dip of 52° and rake of 98°  and extend the fault from the surface to 20 km depth (Arrowsmith et al., 2017).

Almaty City Splay
Using the profiles from satellite-derived DEMs and our field GNSS measurements (Section 2.1), we constrain location and lateral extent of the city splay initially identified by Grützner et al. (2017). This shows that an active fault is cutting through the alluvial fan Almaty is built on; the splay steps into the basin instead of following the edge of the mountain range. We assume that this splay could continue through the city and join with the mountain-range faults. We use a 25 km length fault, consistent with a magnitude 6.5 earthquake (from a Wells & Coppersmith, 1994 relationship), and assume the fault would extend from the surface to 20 km depth. We model it as a purely thrust earthquake (rake = 90°), with a dip of 45° south (lacking any other dip constraints). This fault extends under Almaty districts of Nauryzbay, Bostandyk, and Medeu.

Northern Splay
In Section 2.1, we mapped a fold north of Almaty and conclude this is caused by a thrust fault, dipping south under the city. We model this fault as a 25 km length fault, equivalent to a moment magnitude 6.5 earthquake, dipping at 45° from the surface down to 20 km depth with a rake of 90°. Though we note, due to the narrow width of the anticline the fault may flatten at relatively shallow depth (Section 5.5). This fault extends at depth under many of Almaty's districts (Alatau, Jetysu, Almaty, Auezov, and Nauryzbay). Grützner et al. (2017) mapped folding in Talgar, along the range front to the east of Almaty, and interpreted it as the surface expression of a blind thrust fault. We use a 25 km length fault along Grützner et al. (2017)'s mapped ruptures as an input for an earthquake scenario, with a magnitude 6.5 earthquake. We use a rake of 90°, dip of 45° south with the fault extending from the surface to a depth of 20 km, for consistency with the other scenarios.
10.1029/2021EA001664 11 of 25 PGA for the Chon Kemin earthquake. For their Chilik scenario, Almaty experienced 0.1-0.2 E g PGA, which is higher than the 0-0.1 E g PGA in our scenario. This is due to Mosca et al. (2019) extending their fault much closer to Almaty. Some differences arise due to the different GMPEs used (see Section 2.2).

Damage and Loss
The OpenQuake Engine calculates the losses (fatalities and financial replacement cost) and damage (buildings with no damage, slight, moderate, extensive, and complete) at each exposure point for each different building class. All calculations are performed at "night", that is, assuming that all occupants of a building are home. In order to display this information, we take the average of the three GMPE results at each exposure point. The only exception is in the northern splay scenario, in which we only display the Abrahamson et al. (2014) result as we suggest this shows the most plausible result for this fault geometry (see Section 5.5). We then aggregate the damage and losses for each district in Almaty.
The totals for each scenario are shown in Table 2, and the loss and damage per district is shown in Figure 7. Note that for damage estimates we display the number of completely damaged buildings (i.e., those structurally damaged to such a degree that whilst the building may still be standing, it must be removed and rebuilt). Complete damage does not necessarily mean complete collapse-from a risk perspective this would mean a complete economic loss as the building must be rebuilt, but as a rule of thumb 20%-40% of completely damaged buildings collapse. Because of this and the variety of ways in which a building can collapse, neither a completely damaged building nor a completely collapsed building will mean that all occupants of that building will die, instead these losses are typically in the range of 10%-40% of the occupants. Hence in the following figures, fatality percentages are much lower than collapse percentages.
The Chilik and Kazakh platform scenarios cause the least damage, cost and fatalities, despite being amongst the two largest magnitude scenarios that we ran, primarily due to their large distance from Almaty. The Chon-Kemin also is one of the least damaging of these scenarios, but we estimate it would still cause ∼1,700 ± 2,100 fatalities, 2,500 ± 1,900 million USD structural cost and completely damage 6,500 ± 3,900 buildings if this earthquake were to occur on the current population.
The Northern splay, city splay and 1887 Verny scenarios show the highest damage, cost, and fatalities. This is due to their proximity to Almaty, with northern and city splay faults extending underneath Almaty city itself, causing severe estimated damage. We note that whilst we used a higher magnitude for the Verny scenario, 7.3 compared to 6.5 for the two splays, the damage and losses are lower for Verny, as the magnitude of shaking is highest on the surface projection of the fault (Figure S4), and the Verny fault is under a smaller area of the city.
The district with the most completely damaged buildings is strongly affected by the location of the fault: scenarios to the south of Almaty such as 1887 Verny and 1989 Chon Kemin most strongly affect southern Note. Out of these scenarios, the hypothetical northern splay earthquake would cause the largest loss and damage to the city of Almaty. The uncertainties from the range of GMPEs used are very large in all cases. Abbreviation: GMPE, ground motion prediction equation.

Building Height Results
The risk calculations we perform rely upon exposure models accurately representing the buildings. But this is a constantly evolving feature, as cities expand and grow through time. Here, we compare building heights derived from high-resolution DEMs with field measurements to see if DEMs can reliably estimate city building heights, and compare the relative accuracy of estimates derived from SPOT and Pleiades imagery, with a view to incorporating these into exposure models to help improve risk calculations in the future.
We used Google Earth to identify flat-roofed buildings, to ensure that the height of building we measure in the field matched the potential height a satellite DEM might retrieve ( Figure S7, Table S6). We used a laser range finder (TruPulse 360R) to calculate the vertical height distance between 15 buildings and our location on the ground, measured using GNSS. We measured at least two corners of each building (depending on visibility from the pavement) and repeated each measurement; the "measured" errors in Figure 8 were calculated using the standard deviation of our field measurements.
Some locations at which we made field laser measurements were under tall trees, meaning that the satellite estimate gave implausibly low building heights: overestimating measurement location height leads to underestimating building height when subtracting absolute building height from absolute measurement location height. In light of this, for those locations we chose another, close location on the same road and instead measured the relative building height to this location. For buildings 1 and 2, we chose a location ∼30 m south east, and for building 3 we chose a location ∼55 m south.
We then filtered the satellite-derived pointclouds using polygons of each building from OpenStreetMap. To make our measurements directly comparable, we subtract the satellite-derived height of the building from the satellite-derived height of our location on the street. We used a 5 m buffer around our location and average the points in the pointcloud. The satellite measurement errors in Figure 8 were calculated using the standard deviation of the building points and the measurement location points.
All of the Pleiades scenes were processed using Agisoft Metashape, as we found that no pointclouds created in ERDAS Imagine 2020 recovered the heights of the ∼100 m high buildings 12 and 13 and instead estimated heights around 20 m. We attribute this to these tall buildings erroneously being removed as noise, despite varying processing parameters to attempt to mitigate this. For the SPOT DEM we processed it using ERDAS Imagine, because these tall buildings were not recovered in either software, but the ERDAS Imagine  Markers are slightly horizontally offset from the field-measured value for clarity. Note that the Pleiades west only covers the area in which buildings 11-15 reside. Plotted error bars are one standard deviation. Note that field measured errors are  E 0.6 m and are consequently mostly too small to see. pointcloud produced heights much closer to the field measured heights. We attribute this to a lack of tiepoints between the SPOT scenes found in the city in Agisoft.
No Pleiades Stereo Scene west value is given for building 15, as it is the shadow of buildings 12 and 13. Additionally, the SPOT measurement of building 8 is likely to be aliasing the height of neighboring building 1.
Most (80%-100%) of the separate satellite-derived building heights underestimate the field-measured height ( Table 3). The Pleiades tri-stereo pointcloud performs best, with the lowest average percentage difference (13.6%) and the largest number of buildings measured within error (80%). This result is consistent with expectations (Zhou et al., 2015), as Pleiades has a higher resolution than SPOT and tri-stereo has an additional nadir image to help constrain the retrieval of the building height. SPOT percentage differences and errors are highly influenced by the poorly resolved buildings 8, 12, and 13; excluding these buildings SPOT has an average percentage difference of 16.5% and 75% of buildings measured within error.
We note that this analysis is performed on a small number of building heights and statistics drawn from these measurements should be considered accordingly. But here we show a proof of concept that satellite-derived pointclouds can resolve building heights in Almaty city (discussed further in Section 5.4).

Potential Under-Estimation of Losses
As commented in Section 2.2, the exposure model that we use is potentially biased toward high-rise buildings. Table S2 shows that the total buildings in this exposure model are ∼50,000 for ∼1.6 million occupants. This equates to slightly over 30 people per building, which is a reasonable occupation for high-rise buildings but may under estimate the significant number of 1-2 storey masonry buildings in Almaty. Due to the longer natural periods of vibration, high-rise buildings are usually less prone to earthquake damage (Silva et al., 2015). This may mean that the results we present here are slightly underestimated.

Comparison to Historical Results
We compare the loss and damage from the three historical earthquakes at the time they occurred to the estimated damage if these events were to happen today. Reports from Hay (1888)  Note. Pleiades tri-stereo performs the best, with lowest average percentage difference between field measured results and the pointcloud and the most buildings height estimated within error. SPOT percentage differences and errors are highly influenced by the poorly resolved buildings 8, 12, and 13, without these buildings the average percentage difference is 16.5% with 75% of buildings measured within error. a We stress that Pleiades stereo west only covers five buildings, whereas the others cover 14 or 15, thus making comparison difficult. There was not a specific expedition to study the effects of the 1889 Chilik earthquake , which instead was studied by questionnaire. As such it is less clear the extent of damage and loss for this event, though Krüger et al. (2015) state that it caused fewer fatalities and damage than either Verny of Chilik and that "the damage was not strong; in a few houses the chimneys collapsed." This is consistent with the findings in this study, in which we estimate if Chilik were to occur today there would be 500  E 1000 fatalities (0.03%) and 2,100  E 2,000 completely damaged buildings (4.2%); though these estimates are solely estimates in Almaty city and towns and villages closer to the earthquake would likely experience higher damage and losses. The large errors on this event are due to the high magnitude (M w E 8.0) and large distance of Almaty from the fault ( E 25 km to the nearest fault plane).
In 1911, the Chon Kemin earthquake killed 452 and injured 740 (Delvaux et al., 2001). It is not clear what the population and building stock of Almaty was at this time, after the damage of the Verny earthquake 24 years earlier. We calculate that with the current exposure of Almaty, a Chon-Kemin event would cause 1,700  E 2,100 fatalities (0.1% of the population) and 6,500  E 3,900 completely damaged residential buildings (12.9%). We note these OpenQuake risk estimations do not include any effects of liquefaction, which was reported to have been induced at depth due to a shallow frozen soil layer (Alshembari et al., 2019). Modeling the sensitivity of the season by adjusting Vs30 to compare results from dry soil (summer) and water-saturated soil (spring and water) would be an interesting area of future study.
From the historical records, the Chon Kemin event was the most destructive, followed by Verny and finally Chilik, though it is not clear if all the reported damage and fatalities from Chon Kemin were sustained in Almaty. Here we calculate specifically for Almaty and find that Verny would be the most destructive, followed by Chon Kemin and then the Chilik event.
For all these events, there are also details of numerous landslides, particularly in the mountainous regions, and some of which were local to Almaty (in the example of the Verny earthquake). Any losses from landslide events have not been taken into account in the loss and damage totals.

District and Building Typology Most at Risk
Here we discuss whether any districts are more at risk from earthquakes than others. Districts with more buildings naturally have a larger absolute amount of damage and absolute number of losses for example, Medeu and Alatau districts (Figure 7). In terms of percentage damage, we find that districts perform similarly. For example, in the Chilik and Kazakh platform scenarios in which the fault does not intersect any of the city districts, there is a fairly consistent percentage damage and loss across the districts, slightly increased in the northern districts, likely due to the lower values of Vs30 ( Figure S2). If a district had consistently higher percentage damage and loss than others, it could indicate poorer building materials in that district, but we do not observe this in Almaty. This is consistent with the similar distribution of residential building types across districts ( Figure 5).
Scenarios in which the fault intersects city districts show increased loss and damage in those districts, as expected. The district Bostandyk (south-central) consistently shows high if not the highest absolute number and percentage of completely damaged buildings, economic losses and fatalities. We ascribe this to Bostandyk being directly above the fault plane in many scenarios: Verny, City splay and Northern splay, and it is the second closest district for the Chon-Kemin event.
We conclude in these scenarios that the main influencing factor in a district's risk is proximity to fault: for distal earthquakes the districts behave similarly to each other due to their similar breakdown residential building types. It is primarily which districts are near to or above a fault plane that has the most affect.
AMEY ET AL.

Utility of Calculating Building Heights
In Section 4 we showed that some satellite derived DEMs can accurately capture field-measured building heights. Pleiades tri-stereo performs particularly well with 13.6% average percentage difference between measured and point cloud-derived measurements, and 80% of buildings measured accurately within error.
SPOT is a lower resolution data set than Pleiades, so a larger area can be covered with lower cost. Without the two tallest buildings (12 and 13), SPOT performs comparatively well with Pleiades stereo and tri-stereo, making SPOT the most cost effective data set for this calculation. However, given that buildings 12 and 13 were not resolved at all (when processed in both Agisoft and ERDAS) in SPOT DEMs, in order to have confidence in identifying buildings heights across an entire city and not missing high-rise buildings, likely to house many more inhabitants than smaller residential buildings, we recommend Pleiades tri-stereo as the best data set for calculating building heights.
We noted in Section 4 that the satellite-derived building heights are systemically underestimates of the field-measured building height. This may be in part due to the building in the pointcloud being misaligned to the Open Street Map polygons of building footprint. However, we feel this underestimation is also largely due to the DEM processing. This could be due to filtering of higher elevation points as noise for the buildings, or could be an incorrectly high elevation of the ground due to interpolation between buildings. To improve this, the size of the building polygon could be reduced to limit the number of "ground" points being incorrectly incorporated into the building polygon. Additionally, here we have also simply presented the average of all the points of the pointcloud within a building polygon and one sigma uncertainty, however this could instead be done by more robust method of plane fitting using confidence estimates on each point, or by cluster analysis on the highest points within the building polygon. It should be noted as well that there may be small errors in our field measured estimates.
For cities such as Almaty that have freely available shapefiles of all the buildings (available from: http:// download.geofabrik.de/asia/kazakhstan.html, accessed 8/11/2019) then automatically calculating building heights over the whole city would be an achievable task. This could be achieved by using the cadastral map (map of property information), taking the midpoint of each polygon and assigning to corresponding pixel on a DEM, or to nearest point in pointcloud. If details of the building outlines were not available, machine learning methods present a very viable method of identifying building footprints across an entire city, potentially using a combination of optical imagery and DEMs. Studies have calculated building height by extracting building footprints automatically from Pleiades DEMs and in addition, reconstruct the roof geometry in 3D (Lafarge et al., 2006(Lafarge et al., , 2008, which currently is very time consuming (Lafarge et al., 2008 note computing time is 3 h for a few tens of buildings) but could potentially eventually be applied at city-scale.
Here, we calculate the difference between the elevation of a building (calculated by averaging the pointcloud within a building polygon) and the elevation of the point on the ground where we performed the calculation, in order to be directly comparable to our field measurements. But the building heights could instead be calculated across a city by calculating the difference between a DEM and the DSM derived from satellite imagery, for each building polygon. Additionally, height could be estimated by calculating the difference between a high resolution DEM like Pleiades and one which only captures longer wavelength topography, like SRTM at 90 m resolution.
One difficulty would include buildings whose roofs are not flat. In this study, in order to accurately measure buildings from the pavement we targeted flat-roofed buildings, and thus all DEM measurement points on the building should all be measuring the same height. If instead a building had a pointed roof, or features on top such as domes or multiple levels then the height found would be an average of these. High resolution images such as Pleiades and techniques such as those used by (Lafarge et al., 2008) would likely to be essential to accurately map at large scale.
For the assessment of damage or losses using the OpenQuake Engine, the number of storeys is used instead of building height. Assuming a linear scaling factor would be used to scale between building height and number of storeys, an error in building height would scale linearly with error in number of storeys. An error of  E 3 m that we find here would result in an error of   E 1 storey. Further analysis is required to quantify the affect this would have on risk calculations. AMEY ET AL.

10.1029/2021EA001664
17 of 25 An additional difficulty lies in the automatic calculation of the height of very tall buildings. We found that buildings 12 and 13, which are over 100 m tall, were not resolved at all by ERDAS Imagine 2016's photogrammetry tool. In Agisoft we were able to find the heights for them, but their footprints were slightly misaligned and aliased onto that of the smaller (∼30 m) building at their base. Thus for high buildings, these need to be manually checked.
The advantage of using satellite data to estimate exposure is that it has the potential be updated rapidly on a global scale, without the need for expensive and time-consuming field campaigns (see Pittore et al., 2017 for further discussion). Satellite data is already used to group areas of similar building classes (e.g., Wieland et al., 2012), so the proof-of-concept we present here is a natural extension to exposure classification.

Fault Geometry at Depth
The most damaging scenario modeled here is a M w E 6.5 earthquake on the Northern splay. In Section 2.3.2.2 we parameterized it as a fault that dips uniformly at 45°, but due to the narrow width of the anticline (Figure 4) it is highly likely that this fault may initially dip at ∼45° then flatten at a shallow depth, such as seen in Tabas in Iran (Walker et al., 2003). To demonstrate the importance of properly parameterizing this fault at depth, we compare a fault dipping at 45° and then shallows at 5 km depth to the result of a planar fault dipping solely at 45° to depth.
We find that a fault whose dip shallows out at 5 km depth has 40% more damaged residential buildings, 38% more fatalities in residential buildings, and the cost of replacing the residential buildings is 27% higher (Figure 9). This is a significant increase on what was already the worst scenario modeled here and has large implications for Almaty. There are also a number of other cities in the northern Tien Shan where cities are located near the mountain range front and are potentially at danger of faults stepping out into the basin, underneath the cities, including Bishkek, the capital of Kyrgyzstan (population estimated ∼1 million in 2020 [United Nations, 2018]) and Taraz in Kazakhstan (population estimated 400,000 in 2020 [United Nations, 2018]). As the capital, Bishkek is the economic center of Kyrgyzstan and its proximity to the Issyk-Ata fault potentially poses a significant risk (K. Abdrakhmatov et al., 2003;Bindi et al., 2012;Erdik et al., 2005), and there may be other unidentified faults in the city. More work is required to characterize the earthquake risk to these cities as they continue to grow into major metropolises.
This result particularly emphasizes the need for proper characterization of fault location and geometry, particularly in the Tien Shan, and other areas at risk of buried faults close to major cities. This motivates further studies such as seismic reflection surveys to characterize faults at depth (e.g., Pratt et al., 2002) and investment in seismic networks across regions like central Asia for areas where they do not already exist.
The three GMPEs that we use here calculate very different ground motion fields for a ramp-flat fault geometry due to different hanging-wall terms in their equations ( Figure S5). The Campbell and Bozorgnia (2014) GMPE finds the highest magnitude of shaking at the deepest part of the fault, which would have disastrous consequences as it puts an extremely high shaking extending under much of Almaty. The Abrahamson et al. (2014) and Chiou and Youngs (2014) GMPEs predict almost equal magnitude shaking at the surface portion of the fault as at the deepest section at the fault. In Figure 9 we used Abrahamson et al. (2014)'s GMPE as we suggest this shows the most realistic ground motion field, with a maximum at the surface and a minimum along the fault hingeline. How the GMPEs calculate ground motion fields from geometries at depth has significant implications for characterization of hazard in shortening regions, as many active, reverse faults may show ramp-flat-ramp geometry. This is discussed further in the next section.

Hanging Wall Effects
One difficulty of performing hazard calculations with faults located very close to cities is that most GMPEs are empirical, regression fits to strong-motion data. To accurately represent ground motion there need to have been recordings of ground motion for a given type of earthquake at a given distance, which means there are fewer measurements close the fault. But the ground motions in a fault hanging wall are significantly amplified in a reverse faulting events; short period ground motions may be doubled on the hanging wall, compared to a similar distance on the footwall (Abrahamson & Donahue, 2013). Damage caused by Figure 9. Comparison of hazard and risk for the northern splay, parameterized as dipping at 45° and then shallowing out to a 5° ramp (left column, using the GMPE of Abrahamson et al., 2014), and a fault that dips at 45° (right column). Top of fault shown in red, fault at depth shown in blue, and hatched line indicates hinge-line. Peak ground acceleration, fatalities, damage, and economic cost is at least 25% higher for the ramp-flat fault parameterization. AMEY ET AL.

10.1029/2021EA001664
19 of 25 reverse faulting earthquakes is also more extreme and more spatially variable in the hanging wall than in the footwall (King et al., 2018).
Due to the large number of events and measurements taken in California, where strike-slip faulting on approximately vertical faults is prominent, many GMPEs do not include hanging wall effects. Even GMPEs that do include hanging wall effects, such as those used here (Section 2.2) show a large range in ground motions for the same input fault ( Figure S5). This is because there is insufficient near-fault strong-motion recordings data from reverse faults to properly constrain the empirical GMPE equations. The 1994 Northridge California and 1999 Chi-Chi events are some of the few examples of recordings in both the footwall and hanging-wall (Abrahamson & Donahue, 2013). Additionally, amplification due to directivity varies significantly depending on assumptions into directivity models for reverse faults: the control of the hanging-wall scaling assumptions on directivity is greater than that of the data (Spudich et al., 2014). This is a particular problem when trying to quantify the risk due to faults underneath cities, where hanging wall effects are felt under the city itself, as in our city splay and northern splay examples. In order to better quantify the hazard and risk of cities such as Almaty we need to better constrain the relative hanging-wall and footwall effects. In order to do this, we either need to improve GMPEs with better measurements of earthquake ground motions on the hanging-wall, better understanding of the fault geometries at depth and equations that match these observations as well as having reasonable built-in assumptions, or to instead use a full wave propagation modeling approach. Seismic-wave propagation is more computationally expensive than using GMPEs and relies on accurate velocity structure models. But with good knowledge of parameters, through seismic-wave propagation models the shaking in foreland basins can be better quantified, with parameters such as source depth, fault geometry and basin structure heavily influencing the ground motions (O'Kane & Copley, 2020).

Fault Parameterization-Chon Kemin Example
As shown in Section 5.5, the parameterization of a fault can have significant effect on the expected hazard and loss. To further test the sensitivity of our models, we test different fault geometries for the historical 1911 M w E 8.0 Chon-Kemin earthquake. This rupture has been well mapped by Arrowsmith et al. (2017) and so we are able to parameterize the fault at different levels of complexity, using the mapped surface expression of faulting as a guide. We model this as one-fault (striking 264° and dipping at 52° determined by Kulikova and Krüger (2015)), two faults (fitting the mapped surface ruptures [Arrowsmith et al., 2017] and dipping at 52°) and a complicated, multi-segmented fault using Arrowsmith et al. (2017)'s complex fault network solution ( Figure 10). For each of these we use a magnitude of 8.0, rake of 98°  and assume the fault uniformly dips from the surface to 20 km depth.
The risk results for these three parameterisations ( Figure 11) show that whether the fault is parameterized as one, two or a complex fault rupture, the losses at district level show the same relative rankings: we find that Alatau suffered the most fatalities, followed by Bostandyk district and so on, and this is the same order for each fault parameterization. Whilst the relative losses and damage are consistent between model runs, the two-fault parameterization has consistently higher magnitude of losses and damage than the other fault set-ups. We attribute this to the down-dip extent of the fault plane being nearer to Almaty in the two fault set-up: Figure 10 shows that the two-fault scenario overlaps with Almaty's southern city district (Medeu) whereas the other two fault set-ups do not. A smaller distance between the fault and the city means that Almaty experiences much higher estimated shaking, because the GMPEs use distance to the fault in their calculations and consequently Almaty experiences more loss and damage. We note that the largest source of uncertainty for these calculations comes from the GMPEs and the current limits they have in computing ground-motion in presence of multi-segmented ruptures, and when comparing the damage and losses for each district of Almaty the results of these scenarios are within error of each other ( Figure S6). We also note that many empirical GMPEs are not well-tested in the case of multi-segmented rupture.
The consistently increased damage and loss for a two-fault model shows the importance of properly parameterizing fault geometry. This suggests that in terms of active tectonics research helping seismic hazard assessment, mapping the location of surface ruptures and being able to give estimates of a fault dip, extent AMEY ET AL.
10.1029/2021EA001664 20 of 25 of down-dip rupture and the geometry at depth (as in Section 5.5) are the most useful pieces of information we can provide.

Verny Earthquake With Simulated Historic Exposure
No surface rupture from the 1887 Verny earthquake was described at the time or has since been identified, so it remains unknown how close to Almaty the earthquake occurred. The report of Hay (1888) provides great detail about the effects 1887 Verny earthquake, and so we use this old account of Verny to create a simple exposure model, to test whether the fault model used in Section 2.3.1 for the Verny earthquake gives similar results to the contemporary account (e.g., Villar-Vega & Silva, 2017, in which the authors replicated past events).
To match the account of Hay (1888), we create a simple exposure model comprised of 1,839 masonry unreinforced (i.e., brick) buildings, with 100 fired clay bricks and 1,739 adobe bricks, and additionally 840 wooden houses. We evenly distribute the reported 25,000 inhabitants of Almaty at this time between these 2,679 residential houses. We locate all these building classes at one exposure point at the oldest part of the city, to best represent the location of the population at the time: 76.945 °E, 43.26 °N, in the present-day district of Almaly, less than 100 m west of the district Medeu.
From the Hay (1888) report, the Verny earthquake damaged all the stone buildings and killed 180 from Verny town itself, 328 people in total including the surrounding regions.  (1888).
To test whether our estimates are nearer the reported values, we instead model the Verny earthquake fault as rupturing along the mountain range front, much closer to the city of Almaty than our modeled fault in Figure 10. Three different fault parameterisations (red = fault surface, blue = fault at depth) using the surface Chon Kemin ruptures (orange) mapped by Arrowsmith et al. (2017). E 427 undamaged buildings, which is closer to the Hay (1888) report as again, within error, all the buildings may have been damaged. Based on the building damage, we conclude that both the scenario described in Section 2.3.1, and a scenario rupturing the range front could potentially damage all buildings within error.
The fatalities that we model are significantly lower than those reported at the time: our models estimate tens of deaths, whereas the figure recorded by Hay (1888) was 180 in Verny. The main reason for this underestimation is likely due to the fragility and vulnerability curves that we use. These curves have been developed for modern building classes and hence using these on historical building is very likely to assume higher resilience to earthquake ground motions than is correct for these buildings. This means that the building damaged and fatalities will be underestimated for this historical event. We also lack details from the Hay (1888) report on the causes of these fatalities. The OpenQuake calculations here are only for fatalities caused by the earthquake itself, and the Hay (1888) report total likely includes fatalities due to landslides Figure 11. Left column shows buildings with complete damage, central column shows fatalities and third column shows economic loss for three Chon Kemin scenarios, where row one is a one-fault set-up, row two is a two-fault set-up and the third row is a complex-fault set up (see Figure 10). Whilst the three scenarios show the same relative damage or loss between districts, the two-fault set up has consistently higher values than the other two-fault parameterisations. Blue line indicates the northernmost extent of the buried fault planes. Note the scales are different to Figure 7. ("earthslips" in the report). The report also does not clarify whether all deaths were from the earthquake itself-the reports were from June 10, 1887, almost 2 weeks after the earthquake on May 28th, and it is possible the figures quoted could include people who survived the earthquake itself, but died due to other causes.

Conclusions
In this study, we mapped geomorphic expressions of fault splays from satellite DEMs and ran a suite of earthquake scenarios for the city of Almaty, Kazakhstan, in Central Asia. Almaty has been exposed to several large ( M w E 7.3) earthquakes since it was founded in 1854, but the city has expanded dramatically since these previous events, thus the exposure (and potentially risk) is much greater. Using high-resolution optical stereo satellite imagery (SPOT and Pleiades), we create high resolution DEMs with which we have identified subtle fault splays within Almaty city and under folds in the north. We assess the plausibility of using these satellite-derived DEMs to estimate building heights, with an aim to update exposure models as cities expand. We find the satellite scenes are able to resolve the height of buildings within error satisfactorily (Pleiades tri-stereo estimates 80% of building heights within one sigma uncertainty, with an average percentage difference of 14%). We use the software OpenQuake Engine to assess the damage, fatalities and economic losses for a series of earthquake scenarios. These include three historical earthquakes known to have damaged Almaty, and four plausible earthquakes that we base on our observations of city splays and palaeoseismic evidence from the literature. These scenarios show that the districts with the greatest number of residential buildings, Medeu and Alatau, frequently have highest collapses but the particular district at most risk is highly affected by the exact location of the fault. In our scenarios, the worst case explored would be if a fault to the north of Almaty that dips south under the city were to rupture in a magnitude 6.5 earthquake, causing a total of ∼12,300  E 5,000 building collapses, ∼4,700  E 2,700 Million US dollar cost to be incurred and ∼4,100  E 3,500 fatalities. We also note that the values calculated here may be under-estimates, due to the exposure model used. In all scenarios, districts directly above a buried fault show the highest percentage damage and losses. This, along with our discussion on how fault location and geometry affect the hazard and risk calculations, demonstrate that it is a priority to properly map and characterize the location and geometry at depth of active faults within cities to adequately characterize the risk. This could include through field mapping, seismic reflection surveys and seismic networks.

Data Availability Statement
Pleiades images made available by CNES in the framework of the CEOS Working Group for Disasters. ©CNES (2012©CNES ( , 2014©CNES ( , 2017, and Airbus DS, all rights reserved. Commercial uses forbidden. TanDEM-X mission topography DEM data made available through the German Aerospace Centre (DLR) science service (proposal number DEM_GEOL1336-Active Tectonics of the Northern Tien Shan). ASTER GDEM is a product of METI and NASA. AW3D original data are provided by JAXA. The DEMs derived in this study are available on OpenTopography (datasets https://doi.org/10.5069/G9H41PMP, https://doi.org/10.5069/ G9MS3QZT and https://doi.org/10.5069/G9RJ4GN3). The OpenQuake engine (version 3.8.1 used here) is an open-source, freely available code, see https://www.globalquakemodel.org/openquake to download the latest version. The input datafiles used to run the OpenQuake models along with scripts are available at https://github.com/ruthamey/almaty_hazard_risk.