Scenario ensemble modelling of possible future earthquake impacts in Bhutan

Recent large earthquakes in the Himalaya have resulted in tens of thousands of fatalities, yet these events are thought to have had relatively moderate magnitudes for the region. Evidence suggests multiple events throughout the Himalaya in the last 1000 years have had M > 8.0 and at least two have had M > 8.5. Despite this, understanding of earthquake risk in the region is poorly constrained, particularly in Bhutan, where research on both past and future earthquakes is notably scarce. While recent work has clearly shown the potential for large earthquakes here, the impacts from potential future earthquakes in Bhutan are entirely unknown. This study attempts to address this by modelling the potential fatalities associated with a range of plausible earthquakes through a scenario ensemble analysis in order to inform contingency planning and preparedness. The results show that both the timing and location of future earthquakes are critical factors in determining the number of fatalities, with night-time earthquakes, and those located in the west of the country proving most fatal. The worst case involves ~ 9000 fatalities and results from an M8.5 earthquake directly beneath Bhutan. Nevertheless, at the local scale the number of fatalities appears to saturate at ~ M7.5, since larger earthquakes do not result in significantly larger modelled fatalities. This suggests that local-scale impacts approaching the worst case may be relatively common and emergency planning could focus on comparatively moderate-sized earthquakes since larger, less frequent events may not necessarily result in significantly more fatalities.


Introduction
Earthquakes are amongst the most deadly natural hazards worldwide (Bilham 2004), and populations living in the Himalayan region are amongst the most at-risk to earthquake disasters globally (Robinson et al. 2019). In the last one thousand years, at least 15 large earthquakes are thought to have struck the region (Fig. 1), several of which are likely to 1 3 have had M > 8.0 (Feldl and Bilham 2006;Kumar et al. 2006;Mugnier et al. 2013;Sapkota et al. 2016;Bilham 2019). Due to the large present-day population living in the Himalayan region, recent earthquakes have proven especially deadly, with the 2005 Kashmir earthquake in Pakistan resulting in an estimated 73,000 fatalities (Bilham and Hough 2006) and the 2015 Gorkha earthquake in Nepal resulting in ~ 9000 fatalities (Goda et al. 2015). However, both these earthquakes had relatively moderate magnitudes for the region (M7. 6 and M7.8, respectively) and current estimates suggest that an earthquake with M > 8.0 could result in excess of 150,000 fatalities (Wyss 2005;Robinson et al. 2018). Nevertheless, both the hazard and risk from earthquakes in this region are poorly understood compared to other seismically active regions globally. This is especially true for Bhutan, where research on earthquake hazard and risk is notably sparse. Until recently, Bhutan was the only section of the Himalaya where evidence for a past major earthquake was lacking. Combined with present-day earthquake rates being significantly lower in Bhutan compared to the rest of the Himalaya (Drukpa et al. 2006;Gahalaut et al. 2011), it had been thought that the region was aseismic and consequently earthquake risk here was thought to be low compared to the surrounding Himalaya. However, it has been suggested that the low modern day earthquake rates may be a temporary result of the 1897 Shillong Plateau earthquake south of Bhutan releasing stress in the region (Gahalaut et al. 2011). Furthermore, recent work in south-west Bhutan has identified evidence for a previous earthquake with M 7. 8-8.3 in ~ 17148-8.3 in ~ (Hetényi et al. 2016 (Fig. 1), which appears to be supported by written reports at the time from southern Tibet Fig. 1 Known and inferred earthquakes with estimated magnitudes on the Main Himalayan Thrust (MHT) and surrounding faults in the last 1000 years compared to modern day urban populations. Inset: slip and earthquake magnitude potential on different segments of the MHT based on the time since last event, highlighting the Central Gap in western Nepal (500 + years) and the Bhutan Gap (300 + years). Adapted from  and Bilham (2019) and eastern Nepal of strong shaking in the eastern Himalaya (Bilham 2019). Consequently, it now seems likely that Bhutan can and does experience large earthquakes like those recently experienced in other parts of the Himalaya.
Despite this emerging evidence for potentially large earthquakes in Bhutan, the impact of such an earthquake remains entirely unknown. The only recent earthquake to have significantly impacted Bhutan was the 2009 M6.1 event in eastern Bhutan that resulted in > 1000 buildings requiring major repairs and 11 fatalities (Kayal et al. 2010). A larger, M6.9 earthquake in Sikkim, India, in 2011 also caused significant damage in western Bhutan, with at least 1000 buildings requiring major repairs, although no fatalities were recorded. Evaluating the potential damage and loss of life from future larger earthquakes occurring beneath Bhutan is therefore essential for understanding and planning the necessary emergency response. However, current approaches for seismic hazard analysis are not well tailored to such contingency planning, particularly in regions where previous earthquake data are sparse (Panza et al. 2011;Wang 2011;England and Jackson 2011). Recent work in Nepal, however, has shown that the use of scenario ensembles to assess the variation in impacts across multiple different earthquake scenarios can provide pertinent information for emergency planning that other approaches cannot (Robinson et al. 2018).
This study therefore focuses on an ensemble analysis of potential future earthquakes affecting Bhutan in order to derive the first-ever estimates of potential earthquake impacts for the country. A suite of plausible future large (M > 7.0) earthquakes and their impacts are modelled using openly available population and building fragility data. The impacts from each scenario are then evaluated, and key summary statistics are derived to evaluate variation across the ensemble. The results are intended to provide first-order estimates of the potential scale of impacts from a future earthquake affecting Bhutan in order to inform government and humanitarian-level contingency planning.

Earthquakes in the Himalaya
The Himalayas are the result of the collision between India and Eurasia, with convergence rates of ~ 18 mm a −1 (Bilham 2019). This collision primarily occurs along the Main Himalayan Thrust (MHT), which runs from Afghanistan to Myanmar and is locked to a depth of ~ 30 km. The MHT is known to exhibit a ramp-like structure, with near-surface dip angles of ~ 30° to a depth of ~ 5 km before transitioning to a subhorizontal fault dipping ~ 10° (Pandey et al. 1995). Previous large earthquakes are known to have occurred on every segment of the MHT within the last 1000 years (Fig. 1). The largest was the 1950 M 8. 6-8.8 Assam earthquake Bilham 2019), although the earthquake in western Nepal in ~ 1505 may have been a similar magnitude (Kumar et al. 2006). Six of these events were recorded in the last century (Mugnier et al. 2013;Bilham 2019), highlighting that earthquake records in the Himalaya are incomplete prior to 1900, suggesting that many more large earthquakes are missing from the records. Nevertheless, there are at least two notable seismic gaps where the last major earthquake was several hundred years ago. An ~ 500 km section in western Nepal/eastern Uttarakhand has the longest gap since the last major earthquake at > 500 years, which equates to a slip deficit of > 9 m and a potential earthquake of M8.5 + (Bilham 2019). The next longest gap is on the ~ 400 km section in Bhutan where the last major earthquake was > 300 yrs ago in ~ 1714; however, while this event certainly ruptured a portion of the MHT in western Bhutan, it remains uncertain if rupture propagated into eastern Bhutan as well (Hetényi et al. 2016). Irrespectively, 1 3 there is a minimum slip deficit on the MHT in western Bhutan of ~ 5 m equating to an M8.0 + earthquake, with the possibility for an M8.5 + event if both the western and eastern Bhutan sections were to rupture simultaneously (Drukpa et al. 2012).

Scenario ensemble analysis
Scenario ensemble analysis is a relatively new approach for seismic risk assessments based on methods used in climate and meteorological modelling that attempts to bridge the gap between probabilistic and deterministic approaches (Robinson et al. 2018). Ensemble analysis takes a scenario-based approach to estimate the impacts from specific earthquake events, but considers a large suite of plausible events to evaluate the frequency and variation in the potential impacts. The aim is to identify a range of geologically plausible earthquakes based on a combination of historical and geological evidence of previous events as well as understanding of local fault behaviour and regional tectonics (Robinson et al. 2018). Because this approach focuses on understanding the potential impacts of such events to inform contingency planning, the probability of each scenario is not considered. Instead, the outcome is intended to evaluate how the impacts vary across the entire ensemble and thus identify whether particular impacts are inevitable irrespective of the earthquake that occurs. Accordingly, this approach is tailored to contingency planning as it allows the distribution and range of impacts to be considered while maintaining specific scenarios. It should not, however, be misconstrued with more widely used probabilistic or deterministic seismic hazard analyses that primarily focus on understanding earthquake hazard rather than impacts.
The ensemble approach is, however, limited in its treatment of the probability associated with each scenario as well as distinguishing between plausible and implausible events. Nevertheless, the approach has recently been successfully applied in Nepal where it is currently informing both the Government of Nepal and the United Nations contingency planning efforts for future earthquakes (Robinson et al. 2018).
This study focuses on an ensemble analysis of potential future earthquakes affecting Bhutan in order to derive the first-ever estimates of earthquake impacts for the country. The results are compiled for each of Bhutan's Dzongkhags (Administrative Division 1), by evaluating the frequency and variation in impacts combined with estimates of worst-case and average impacts across the entire suite of modelled scenarios.

Earthquake scenarios
The first step is to determine an ensemble of plausible earthquake scenarios. The MHT can be split into multiple different segments, largely based on the spatial extent of previous large earthquake ruptures (Figs. 1 and 2). In Bhutan, the MHT may be split into two segments: the western Bhutan segment and the eastern Bhutan segment (Drukpa et al. 2012). The boundary between these two segments is considered to be near the eastern margin of the ~ 1714 rupture (Drukpa et al. 2012;Hetényi et al. 2016), which coincides with a sudden decrease in fault coupling on the MHT in eastern Bhutan (Marechal et al. 2016). This suggests that the 1714 earthquake may have only ruptured the western Bhutan segment, making the eastern Bhutan segment the only segment without a known earthquake (Fig. 1).
Nevertheless, fault motion, dip angles, depth ranges and slip rates on both segments are largely identical (Drukpa et al. 2012) and the decreased coupling in eastern Bhutan may be a short-term phenomenon, possibly linked to the 1897 Shillong Plateau earthquake (Gahalaut et al. 2011). Consequently, both these segments are considered capable of sustaining large earthquakes and the similarities between them means they may be capable of rupturing simultaneously, producing a larger earthquake than would occur if they ruptured independently ( Fig. 1) Drukpa et al. 2012;Bilham 2019).
Despite only the ~ 1714 earthquake being definitively identified on the MHT in Bhutan, the structure and behaviour of the MHT here differ little from elsewhere in the Himalaya (Fig. 2). Consequently, it is considered plausible that similar earthquakes to those recorded elsewhere on the MHT (  While the MHT forms the plate boundary and is therefore the dominant fault in the region, other subsidiary faults are also present and capable of generating large earthquakes. The most notable are the Oldham Fault and the Dauki Fault, situated at the northern and southern margin of the Shillong Plateau, respectively (Fig. 2). These faults are inferred to be steeply dipping conjugate reverse faults that accommodate 'pop-up' of the Shillong Plateau (Bilham and England 2001). The Oldham Fault is considered the most likely candidate for the 1897 M8.2 Shillong Plateau earthquake (England and Bilham 2015), which ruptured from 9 to 30 km depth and is likely close to a maximum event for this fault. While no previous ruptures have been attributed to the Dauki Fault, the obvious expression of the fault at the surface suggests that it is active and capable of sustaining large events (Ferguson et al. 2012), possibly up to M8.0, similar to the 1897 event. While the 1897 earthquake on the Oldham Fault was comparatively recent, a recurrence of a similar event cannot be discounted and therefore this study considers scenarios across the full width of the fault, including the 1897 rupture zone.
Recordings of small (M < 4.0), comparatively deep (30-60 km) earthquakes in the region (Sharma et al. 2018) provide evidence of two further fault zones that do not reach the surface. These small events show lateral fault motion and are located at a greater depth (> 30 km) than the MHT (Kayal et al. 2010;Sutar et al. 2017;Sharma and Baruah 2017;Sharma et al. 2018;Grujic et al. 2018). The western fault zone is known as the Dhubri Fault zone and coincides with the eastern and western extents of the 1934 and 1714 earthquakes on the MHT, respectively. It therefore marks the location between the Everest and western Bhutan segments of the MHT (Fig. 2) and may form an impenetrable barrier that stops earthquakes rupturing on both segments simultaneously (Drukpa et al. 2012;Grujic et al. 2018). Likewise, the Kopili Fault zone in the east marks the boundary between the eastern Bhutan segment and the western Arunachal segment of the MHT where a significant (20°) change in strike occurs, and may also represent an impenetrable barrier to rupture ( Fig. 2) (Drukpa et al. 2012). Seismicity on both these fault zones is recorded between depths of 30 and 60 km, with depths increasing progressively northwards (Sharma and Baruah 2017;Sharma et al. 2018). Both fault zones clearly have potential for large earthquakes, with the Dhubri Fault zone possibly being responsible for an M7.1 earthquake in 1930 (Kayal et al. 2010;Grujic et al. 2018).
Finally, the Yadong-Gulu rift is a normal fault that extends south from central Tibet and follows the western border of Bhutan, where it appears to end near to where it would intersect the Dhubri Fault zone ( Fig. 2) (Drukpa et al. 2012). Historical and instrumental records suggest that this fault has been strongly active along its northern segments, with at least 14 known earthquakes with M > 6.0, the most recent of which was in 2008, ~ 150 km north of Bhutan (Wu et al. 2011). However, understanding of the fault near Bhutan is limited and in this region the fault is likely to be at significant depth compared to its northern segments. Other active faults are known within and near to Bhutan; however, these are not considered capable of sustaining earthquakes with M > 7.0.
Earthquakes are modelled on each of the fault segments shown in Fig. 2 at 0.5 M intervals starting at M7.0 and increasing to each segment's maximum plausible magnitude (Mmax; Fig. 2). Shaking intensities are calculated using the ground motion prediction equations from Abrahamson and Silva (2008). For each segment of the MHT, Oldham and Dauki faults, earthquakes with M7.0 and M7.5 are modelled occurring at both the near surface (< 10 km deep) and at depth (> 15 km deep), while earthquakes with M ≥ 8.0 are assumed to rupture the full width of the corresponding fault. In total, this yields an ensemble of 55 potential large earthquake scenarios for Bhutan (Fig. 3).

Population exposure
Earthquake fatalities are predominantly caused by the collapse of buildings during shaking (Frolova et al. 2011;Jaiswal et al. 2011b;So 2016). In order to estimate the impacts from each of the earthquake scenarios, an assessment of the likely population indoors during shaking is required. While the indoor population is expected to vary considerably at hourly, daily and seasonal timescales, the scale of this variation is currently poorly understood. The limited studies that have investigated the effect of the timing of an earthquake in terms of consequent impacts have, however, shown significant differences between daytime and night-time earthquakes (Scawthorn 2011;Robinson et al. 2018). This effect is expected to be larger in rural regions where most of the population work outside during daylight hours, compared to urban populations who are more likely to work indoors. In countries with mainly rural populations like Bhutan, the difference in total impacts for the same earthquake occurring at night or during the day may therefore be substantial. Consequently, this study models each of the 55 earthquake scenarios for both daytime (midday) and nighttime (midnight) indoor populations.
Population data for Bhutan are taken from the 2017 Population and Housing Census of Bhutan (PHCB) and evaluated at the Chiwog level (Administrative Level 3), the smallest administrative level available. This shows that in 2017, Bhutan had a total population of 721,155 of which ~ 16% (114,551) lived in the capital city Thimphu. With no available data on the population's daytime and night-time locations, expert elicitation was used to generate an estimate of the indoor population on an average day at both midday and midnight. This was undertaken in coordination with a series of experts from the Royal Government of Bhutan and humanitarian sectors during a one-day workshop in Thimphu in June 2019. This concluded that the indoor night-time population was likely to be ~ 100% of the Chiwog's population (Fig. 4). The daytime indoor population, however, likely varies depending on: a. The total number of people employed in the agriculture sector; b. The total number of people aged 65 and over; and c. The total number of children aged 5-14 attending school.
At midday on an average day, all people employed in agricultural work as well as those aged 65 and over are expected to be outside. In addition to this, it was expected that all children aged 5-14 not attending school would also be outside during the day. All other people were considered most likely to be indoors. The PHCB includes information on population age, allowing an accurate evaluation of the population aged > 65 and 5-14 for each Chiwog. However, data on employment type and school attendance rates were only available as percentages of the population at Gewog level (Administrative Division 2). These Gewog-level statistics were therefore used to derive estimates of the population working in agriculture and attending school for each Chiwog to derive estimates of the total indoor daytime population for each Chiwog (Fig. 4). These results show clearly that the differences in night-time and daytime indoor population are most stark in eastern Bhutan.

Building vulnerability
While the majority of earthquake fatalities result from building collapse, some construction materials are more vulnerable to collapse than others (Jaiswal et al. 2011a). Accounting for the different building types across Bhutan is therefore essential. The PHCB contains detailed data on the primary wall construction material for each building, showing six generalized building types across Bhutan ( However, while data on building damage in Bhutan from the 2009 and 2011 earthquakes were collected, the sample sizes are insufficient to generate reliable fragility curves specific to Bhutanese construction styles. Consequently, proxy curves are required that can be considered a best approximation for Bhutanese buildings. While some building types are unique to Bhutan (e.g. rammed earth structures), most reinforced concrete, mud block (adobe), stone and brick structures are similar in construction style to those found in neighbouring Nepal. In the last ~ 100 years, Nepal has experienced at least five damaging earthquakes for which sufficient data on building damage are available to derive local fragility curves (Guragain 2015;Chaulagain et al. 2018;Gautam et al. 2018). Given the similarities in construction styles, this study therefore uses the fragility curves from Nepal as proxies for reinforced concrete, adobe, stone and brick structures in Bhutan (Fig. 6). However, stricter enforcement of building codes for reinforced concrete structures in Bhutan compared to Nepal means that the fragility curve for these structures is likely to be a conservative estimate (Gautam et al. 2016). For wooden and other building types, only fragility curves derived from combined global samples from HAZUS (Jaiswal et al. 2011a) are available. This study therefore uses the curve corresponding to W1 (wood, light frame, < 5000 sq. ft) to represent wooden buildings. For other buildings, which likely span a large range in construction styles, the RM1M (reinforced masonry bearing walls with wood diaphragms) curve is selected, as this represents comparatively moderate building performance (Fig. 6).
These fragility curves represent 'complete damage', which refers to buildings that are damaged beyond repair; however, this does not directly translate to building collapse. Since most earthquake fatalities result from buildings that collapse during shaking, estimates of the percentage of buildings experiencing complete damage that will actually collapse are applied based on global empirical data (Table 1). Similarly, while not all buildings suffering complete damage will collapse, not all people inside collapsed buildings will be killed. Building type appears to play an important role in fatality rates, with lighter materials resulting in significantly fewer fatalities than heavier materials (So 2016). A global Fragility curves for different Bhutanese building types suffering complete damage, where complete damage is defined as a building being beyond repair; however, the structure may not immediately collapse. Curves for reinforced concrete, adobe and stone buildings are taken from Gautam et al. (2018), brick buildings from Guragain (2015), and curves for wooden and other buildings are taken from HAZUS (Jaiswal et al. 2011a) and correspond to building types W1 and RM1M, respectively empirical study by So (2016) identified the range of typical fatality rates for different building types and provided suggested rates to use in earthquake fatality models (Table 1).

National-scale impacts
Out of the 110 scenarios modelled, 65 (~ 60%) result in one or more fatalities somewhere in Bhutan and 43 (~ 39%) result in > 50 fatalities. Of these, 15 (23% of fatal scenarios, 13% of all scenarios) result in more than 1000 fatalities and 5 result in more than 5000 fatalities (Fig. 7). The worst case results from an M8.5 earthquake rupturing both Bhutan segments of the MHT at night, which is modelled to cause ~ 9000 fatalities; however, if this same earthquake occurred during the day, the number of fatalities is modelled to be ~ 5500. Strikingly, an M8.0 earthquake occurring at night on the western Bhutan segment of the MHT is modelled to result in more fatalities (~ 6750) than this daytime M8.5 earthquake, highlighting both the importance of the timing of the earthquake and the level of exposure in western Bhutan. Night-time earthquakes on average are modelled to be approximately twice as fatal as if they occurred during the day (Fig. 7). In general, smaller magnitude earthquakes occurring at night are modelled to result in more fatalities than larger earthquakes in the day. For example, while a daytime M8.0 earthquake in western Bhutan has ~ 4250 modelled fatalities, a night-time M7.5 earthquake in western Bhutan is modelled to cause between ~ 5000 and ~ 5500 fatalities, depending on whether it occurs at depth or at the near surface. Furthermore, there are four different night-time earthquake scenarios that result in > 5000 modelled fatalities, while the only daytime scenario with such numbers has M8.5.
Earthquakes occurring on only the western Bhutan segment of the MHT are significantly more damaging than earthquakes occurring on only the eastern Bhutan segment (Fig. 7). The average number of fatalities from all earthquakes on the western Bhutan segment is ~ 1000 with an inter-quartile range (IQR) of ~ 150-3000. In comparison, the eastern Bhutan segment has an average of ~ 175 and an IQR of ~ 60-1250. This is especially clear when comparing identical earthquakes on both segments. An M8.0 earthquake in eastern Bhutan results in ~ 2500 night-time and ~ 1500 daytime fatalities, compared to ~ 6750 and ~ 4250 fatalities, respectively, for the same earthquake in western Bhutan. Moreover, smaller earthquakes occurring in western Bhutan are modelled to result in more damage than much larger earthquakes in eastern Bhutan; M7.5 earthquakes in western Bhutan produce at least twice as many modelled fatalities as an M8.0 earthquake in eastern Bhutan. This is most likely a direct consequence of the larger population in western Bhutan, particularly in Thimphu Dzongkhag (Fig. 4), but is also likely influenced by the larger number of adobe buildings in western Bhutan compared to eastern Bhutan (Fig. 5) since these are by far the most lethal building type (Table 1). Of the 68 earthquake scenarios modelled occurring on fault segments outside Bhutan, just 7 (10%) result in > 50 modelled fatalities in Bhutan, and all of these are on the Everest or western Arunachal segments of the MHT (Fig. 7). None of these scenarios result in > 1000 modelled fatalities in Bhutan, with the most resulting from an M8.5 earthquake on the Everest segment, which has ~ 750 fatalities at night and ~ 500 during the day. Previous work, however, has suggested that such an earthquake could result in > 120,000 fatalities in Nepal (Robinson et al. 2018). Earthquakes on the Everest segment with M ≥ 8.0 produce more modelled fatalities in Bhutan than any M7.0 earthquake on the eastern Bhutan segment and most of the M7.0 earthquakes on the western Bhutan segment, highlighting that large, distant earthquakes still pose a notable risk to Bhutan, particularly those in the west. Earthquakes on the western Arunachal segment are not as deadly for Bhutan as  Fig. 3 for scenario codes those on the Everest segment, with an M8.5 earthquake in western Arunachal resulting in less fatalities in Bhutan than any M7.5 earthquake on the western Bhutan segment (Fig. 7). The only non-MHT earthquakes that produce fatalities in Bhutan result from earthquakes on the central and northern segments of the Dhubri fault system and the northern segment of the Kopili fault system. However, these result in < 50 modelled fatalities, which is likely smaller than the model's accuracy, and therefore the potential for fatalities from earthquakes on these faults remains uncertain. Earthquakes modelled on the Oldham and Dauki faults and the southernmost segments of the Dhubri and Kopili fault systems appear insufficient to generate fatalities in Bhutan.

Dzongkhag-scale impacts
The number of fatal scenarios varies considerably for each Dzongkhag: Trashigang has the most with 27, while Gasa has the fewest with just 8 (Fig. 8). However, despite having the largest number of fatal scenarios, Trashigang has only five with > 250 fatalities and none with > 500 fatalities. In comparison, Wangdue has 14 scenarios with > 250 fatalities of which 8 have > 500, and Punakha has 12 with > 250 fatalities of which 9 have > 500. In fact, ~ 60% of the earthquake scenarios that cause fatalities in Wangdue result in > 250 fatalities, while 45% of fatal scenarios in Punakha result in > 500 fatalities, suggesting that when impacts do occur in these Dzongkhags, they are usually large. Nevertheless, of the 20 Dzongkhags in Bhutan, only five experience > 500 fatalities in one or more scenarios, while half of all Dzongkhags never experience > 250 fatalities (Fig. 8).
Evaluating the range of fatalities for each Dzongkhag allows the potential variation in future fatalities to be assessed. Since it is impossible to know which of the modelled scenarios is likely to affect Bhutan, it is important to understand how variable the impacts could be in order to inform contingency planning. Thimphu Dzongkhag has the largest range in modelled fatalities, ranging from 1 to ~ 1600; however, Punakha has the largest D a g a n a G a s a P u n a k h a T s ir a n g W a n g d u e B u m t h a n g S a r p a n g T r o n g s a Z h e m g a n g L h u e n t s e M o n g g a r P e m a G a t .
S a m d r u p J o n g .
T r a s h ig a n g T r a s h i Y a n g t s e median fatalities at 440, compared to 133 for Thimphu (Fig. 9). Similarly, Wangdue, Haa and Paro all have larger medians than Thimphu despite having smaller worst cases and narrower ranges. Because Thimphu has a larger number of fatal scenarios (Fig. 8), this suggests that fatalities in Punakha, Wangdue, Paro and Haa are mainly caused by proximal earthquakes, while more distal earthquakes do not result in fatalities. In comparison, while Thimphu experiences large numbers of fatalities in proximal earthquakes, more distal earthquakes are also capable of producing fatalities, albeit in much smaller numbers, likely as a result of the much larger population here compared to other Dzongkhags (Fig. 4). Notably, all Dzongkhags exhibit a positive skew in the distribution of fatalities, indicating that typically fatalities tend towards the minimum (Fig. 9). Chhukha and Samtse have the largest skew values, indicating that the distribution of fatalities in these Dzongkhags is comparatively long-tailed. In comparison, Haa and Punakha have the smallest skew values, with Haa in particular having a close to symmetrical distribution. On average, skew is smallest in Dzongkhags in the West-Central region, where Wangdue also has a markedly small skew value, and highest in the West region, where Thimphu has a particularly notable large skew. Skew values do not appear to correlate with population since Gasa, which has by far the smallest population in Bhutan (3952), has a value comparable to much more heavily populated Dzongkhags.

Number of Fatal Scenarios
An important consideration for earthquake planning in Bhutan is the link between earthquake magnitude and fatalities. Addressing this allows planners to evaluate whether planning should focus only on very large earthquakes, or whether smaller (but still large) earthquakes are also important for planning. Comparing the maximum modelled fatalities in Dzongkhags for each earthquake magnitude shows that in most cases the largest impacts from an M7.5 earthquake are up to 75% as large as each Dzongkhag's worst-case scenario, which always results from an M8.5 earthquake (Fig. 10). In all Dzongkhags in western-central Bhutan (except Gasa), and Haa and Paro in western Bhutan, the largest fatalities from an M7.5 earthquake are > 90% of the absolute worst case from an M8.5 earthquake. This suggests that at a local scale, M7.5 earthquakes can generate fatalities in individual Dzongkhags comparable to much larger earthquakes, and therefore at a national scale, M8.0 and M8. 5  a g a n a G a s a P u n a k h a T s ir a n g W a n g d u e B u m t h a n g S a r p a n g T r o n g s a Z h e m g a n g L h u e n t s e M o n g g a r P e m a G a t .
S a m d r u p J o n g . T r a s h ig a n g T r a s h i Y a n g t s e a g a n a G a s a P u n a k h a T s ir a n g W a n g d u e B u m t h a n g S a r p a n g T r o n g s a Z h e m g a n g L h u e n t s e M o n g g a r P e m a G a t .
S a m d r u p J o n g .
T r a s h ig a n g T r a s h i Y a n g t s e a g a n a G a s a P u n a k h a T s ir a n g W a n g d u e B u m t h a n g S a r p a n g T r o n g s a Z h e m g a n g L h u e n t s e M o n g g a r P e m a G a t .
S a m d r u p J o n g .
T r a s h ig a n g T r a s h i Y a n g t s e

Uncertainties
Earthquake fatality modelling is affected by numerous uncertainties, many of which are unquantifiable and/or irreducible given current knowledge and data availability. For instance, the location and extent of earthquake scenarios are based on known faults and fault behaviour. However, much of the Himalaya lacks good fault mapping, and earthquakes can occur, and have occurred, on buried or previously unknown faults. There may therefore be other earthquake scenarios not considered in this study that could cause largescale fatalities. Further, the scenarios consider only losses from shaking-induced building collapse and do not include potential impacts from secondary hazards such as landslides, flooding or public health hazards which could exacerbate and extend the modelled impacts (e.g. Budimir et al. 2014). As a result, the resolution of the fatality modelling in this study is conservatively assumed to be one order of magnitude, meaning that scenarios or locations with less than ~ 50 fatalities are likely to be unreliable. Scenarios or locations with more than ~ 50 fatalities are more reliable; however, the numbers of fatalities presented in this study should be considered as broad estimates only.

Earthquake impacts in the wider Himalaya
There are at least 43 different earthquake scenarios that would cause fatalities somewhere in Bhutan, of which at least 15 could result in more than 1000 fatalities in Bhutan. The biggest risk exists from earthquakes on the MHT, which is responsible for all of the fatal scenarios with > 50 modelled fatalities. While the results in this study only consider fatalities in Bhutan, the majority of the modelled scenarios are likely to produce fatalities in India, China and Nepal as well. An M8.5 earthquake on the Everest segment of the MHT results is several hundred fatalities in Bhutan (Fig. 7), but previous work has suggested that such an earthquake could result in > 120,000 fatalities in Nepal (Robinson et al. 2018) and possibly similar numbers in Sikkim, India. This means the vast majority of international aid and response to such an earthquake is likely to be focused in eastern Nepal and Sikkim. Likewise, the worst-case scenario for Bhutan involves an M8.5 earthquake on both Bhutan segments of the MHT and could result in ~ 9000 fatalities nationwide. However, previous modelling has suggested that such an event could also cause > 30,000 fatalities in eastern Nepal (Robinson et al. 2018), meaning that impacts would be larger in Nepal despite the earthquake being centred in Bhutan. Very large (M > 8.0) earthquakes on the MHT in this region therefore present a substantial issue for both local and international emergency response since the impacts will be spread across multiple countries. This is particularly problematic for Bhutan, as the small population means impacts in Nepal and India could be substantially larger than those in Bhutan even for earthquakes centred in Bhutan.

Implications for emergency planning in Bhutan
In terms of emergency planning at the local Dzongkhag level, it is apparent that earthquakes with moderate magnitudes (M7.5) can produce similar numbers of fatalities as very large earthquakes. This means that at a national level, earthquakes with M ≥ 8.0 are only worse than earthquakes with M7.5 because they affect multiple Dzongkhags. But for each individual Dzongkhag the impacts from an M8.5 earthquake in Bhutan and a local M7.5 earthquake are broadly equivalent. Because M7.5 earthquakes occur more frequently than M8.5 earthquakes (approximately 1 per 100 years compared to 1 per 500 years in Fig. 11 Cartograms showing the number of modelled fatalities per Gewog (Administrative Division 2) for night-time occurrences of an M8.5 earthquake on both Bhutan segments of the MHT (a) and an M7.5 earthquake on the western Bhutan segment of the MHT (b). Inset: shaking for both events as per Fig. 3 for reference the Himalaya; Feldl and Bilham 2006), local-level planning could focus on preparing for an M7.5 earthquake since this equates to at least 75% of the local impacts from an M8.5 earthquake.
Planning for an M7.5 earthquake may also be suitable at national level. While an M8.5 earthquake on both Bhutan segments of the MHT produces nearly double the national fatalities of the worst M7.5 earthquake, which occurs on the deep section of the western Bhutan segment, the spatial distribution of impacts is broadly similar (Fig. 11). In fact, for western Bhutan, there is little difference in impacts between the two scenarios: fatalities are primarily concentrated in Paro, Thimphu, Punakha and Wangdue, and the total number of fatalities in these Dzongkhags is broadly the same in both scenarios. For instance, in Thimphu Dzongkhag the night-time M8.5 earthquake produces ~ 1500 fatalities, while the nighttime M7.5 earthquake produces ~ 1200 fatalities. The only difference is an M8.5 earthquake also produces fatalities in eastern Bhutan, although the numbers are significantly smaller than those in western Bhutan. At the national scale therefore, planning could focus on an M7.5 earthquake in western Bhutan since these sized events occur more frequently than M8.5 events and, at the national scale, the impacts would be broadly the same should a larger earthquake occur. At local level, Dzongkhag specific plans would account for the lack of fatalities in eastern Bhutan in this national plan, ensuring that this region is not ignored.
While this work has focused on the fatalities resulting from a future earthquake to inform planning, a major factor that has not been considered is post-earthquake logistics. Bhutan is a highly mountainous and remote country, and road access is generally difficult nationwide. Associated with almost all the modelled scenarios in this study will be major landslides (Nowicki Jessee et al. 2020) that are expected to block roads and isolate numerous communities. During the 2015 M7.8 earthquake in Nepal, at least 20,000 landslides were triggered (Williams et al. 2017;Roback et al. 2017) across the affected area and most major roads were damaged and blocked to some degree, severely hampering access to remote regions. This is likely to be similar in Bhutan but with the added issue that getting resources into the country by air is also complicated by the short runway and difficult approach at Paro International Airport (Wangdra et al. 2018). Further studies are required to assess the likely impact of future earthquakes on logistics, but this is expected to prove a major challenge in the response to any large earthquake in Bhutan.

Conclusions
Bhutan has previously been considered to have low earthquake risk compared to the rest of the Himalaya. However, recent work has shown this to be inaccurate, with conclusive evidence showing very large earthquakes can occur and have occurred here in the recent past. The potential impacts of future large earthquakes on Bhutan have, however, remained unknown. This study has sought to address this for the first time by using a scenario ensemble approach to model fatalities associated with possible large earthquakes that could affect Bhutan to inform contingency planning and preparedness. The results show that while Bhutan's relatively small population limits fatalities within the country compared to recent large earthquakes elsewhere in the Himalaya, there remains the potential for up to 9000 or more fatalities in Bhutan. However, both the timing and location of the earthquake are critical factors in controlling the number of fatalities, with night-time earthquakes being twice as fatal as their daytime equivalents. Earthquakes that affect western Bhutan are, on average, three times more fatal than comparable earthquakes in eastern Bhutan. In fact, smaller magnitude earthquakes in western Bhutan are capable of producing more fatalities than much larger earthquakes in eastern Bhutan, highlighting that risk is primarily concentrated in the more heavily populated western areas.
For the majority of Dzongkhags, the largest number of modelled fatalities from an M7.5 earthquake accounts for > 75% of a Dzongkhags absolute worst-case fatalities. Similarly, at the national-scale, an M7.5 earthquake in western Bhutan could result in comparable numbers of fatalities to an M8.5 earthquake affecting the entire country. This suggests that impacts from a future large earthquake in Bhutan could approach the worst case at both local and national scales even in relatively moderate-sized earthquakes. Given M7.5 earthquakes are expected to occur approximately once per century in the Himalaya, planning in Bhutan could focus on such an event given that a larger earthquake may not result in significantly more fatalities. Nevertheless, large earthquakes in Bhutan present a notable challenge for both local and international responses as large impacts are anticipated to occur in multiple countries, particularly Nepal and India. The results of this study combined with recent similar work in Nepal suggest that the number of fatalities there could be up to three times higher than those in Bhutan, even for earthquakes located in Bhutan. Combined with difficult and limited access into the country, large earthquakes in Bhutan are likely to prove particularly challenging for emergency responders.