Probabilistic tsunami hazard analysis for western Makran coasts, south-east Iran

Makran subduction zone, along the southern coasts of Iran and Pakistan, has a wide potential seismogenic zone and may be capable of generating large magnitude (M ~ 9) tsunamigenic earthquakes. Considering ambiguities exist in tsunamigenic source characterization for subduction megathrusts like Makran, where detailed geologic, seismic, and geodetic data are insufficient, the probabilistic tsunami hazard analysis (PTHA) is the most prevalent approach to handle uncertainties and estimate more reliable tsunami hazard. Here, PTHA is performed for the coastal region of the western Makran, south-eastern Iran. Using the logic tree approach, we have considered the uncertainty of maximum seismic magnitude, earthquake occurrence model, continuity of seismic zone, seismic coupling coefficient, rupture depth, presence or absence of splay faults, fault locations, and fault slip distribution in PTHA calculation for the western Makran region. We have derived uniform tsunami hazard maps for two return periods of 475 and 2475 years for two confidence levels, the mean and the 84th percentile. Ground subsidence effects are also evaluated in a probabilistic manner. According to the PTHA results, Chabahar and Sirik towns are at the highest and lowest tsunami risk, respectively.


Introduction
Taking into account the catastrophic consequences of recent large magnitude (M ~ 9) tsunami-generating earthquakes (i.e. the 2004 Indian Sumatra-Andaman and the 2011 Japan Tohoku earthquakes); it is crucial to perform risk analysis studies for tsunami-prone regions. The first step for a risk analysis study is to conduct a tsunami hazard assessment to evaluate coastal tsunami heights from tsunamigenic source models.
There are two main approaches to perform a tsunami hazard assessment: deterministic and probabilistic. The deterministic approach is usually conducted for a formidable scenario that can pose the largest hazard to the region of interest. However, the deterministic technique might cause unanticipated errors, such as determining an inappropriate worstcase scenario (i.e. underestimated or overestimated). On the other hand, from an economic point of view, the assumption of a biggest or worst possible scenario might not be suitable for conventional design goals or strategies, as buildings and facilities have limited lifespans. Thus, a probabilistic tsunami hazard analysis (PTHA) considering in a comprehensive manner both epistemic (associated to our limited knowledge about things which are constant in time) and aleatory or inherent (associated to natural variations that occur between different events) uncertainties (Abrahamson and Bommer 2005) becomes a better option for a tsunami hazard assessment.
The Makran subduction zone (MSZ) (Fig. 1), formed by northward motion of Arabian plate beneath the Eurasia plate, has a historical catalogue of large earthquakes and associated tsunamis (e.g. the Makran tsunami of 1945 with death toll near 4000).
Different seismologic, geologic and tectonic considerations suggest the segmentation of the MSZ into at least two segments, western and eastern. From the view point of seismic activity, almost all historical large to great earthquakes occurred along the eastern part. The eastern segment hosted several historical large magnitude events at 1765, 1810 and 1945 (Ambraseys and Melville 1982), but the western segment has been quiescent at least in the past 300 years. An ambiguous report of a large earthquake at 1483 in the western Makran was criticized by Musson (2009) and attributed to moderate shallow crustal earthquakes in Fig. 1 Location of the MSZ and its seismicity. The pentagons show the historical earthquakes in this region. The circles show the occurred earthquakes from 1900 to 2019. The radius of a circle presents the magnitude of an earthquake. The colour shading indicates the focal depths in km the neighbourhood regions. The eastern Makran is relatively more active than the western one based on the pattern of offshore modern seismicity (Rashidi et al. 2020b). This could be because of the difference in the plate convergence in the west (continental crust collision) and east (oceanic crust subduction) (Mokhtari et al. 2008). Another argument suggesting this segmentation could be the difference between the offsets of the volcanic arc in the east and west (Byrne et al. 1992).
Moreover, it is stated by Rajendran et al. (2013) that "the Sonne fault possibly divides the Makran subduction front into western and eastern segments". There are also some differences between the convergence rate of the western and eastern Makran (Rajendran et al. 2013;Ghadimi et al. 2021, submitted).
Although it is usual to divide the MSZ into two distinct segments, the possibility of a large earthquake on the western segment or a huge event resulting from a cascading rupture across two segments cannot be ruled out (Musson 2009). Penney et al. (2017) showed that the subduction interface in the western Makran may be locked, accumulating elastic strain, and move in megathrust earthquakes (see Frohling and Szeliga 2016a, b;Pajang et al. 2021), though its seismic coupling is reported to be weak (Ghadimi et al. 2021, submitted).
Recently, probabilistic tsunami hazard assessment of the MSZ has been a subject of interest for tsunami researchers and there are several PTHA for the Makran coastlines. A list of these studies and their specifications is summarized in Table 1. They have been conducted at different scales of global, regional and local, and different levels of detail were applied to describe the subduction geometry, fault slip distribution and other technical factors. A major drawback of some of these studies is that the seismicity rate was used instead of long-term fault slip rate to constrain the recurrence rates in the Gutenberg-Richter (GR) model (see Sect. 2.2). Moreover, they ignored the important uncertainties regarding the recurrence model or maximum magnitude and updip rupture depth.
The first PTHA research for Makran was conducted by Burbidge et al. (2009). They produced probabilistic tsunami hazard maps for the Indian Ocean countries, including Iran, considering the effects of tsunamis generated by the Makran, Sumatra-Andaman, and South Sandwich subduction zones. They presented their results as annual probabilities of exceeding certain tsunami wave amplitude in the water depth of 100 m. Although not explicitly mentioned, it seems that the coupling factor in this study has been taken into account with the technique in which the rate of earthquake occurrence for each subduction zone is determined based on the percentage of the global subduction seismicity. Activity rates for each subduction zone, defined as the fraction of the global subduction zone seismicity that would be expected on each zone, were calculated based on the ratio of its length to the length of all tsunamigenic subduction zones in the world. The main problem in this technique is the assumption that the coupling factor for any given subduction zone is not significantly different from the global average along megathrusts In Burbidge et al. (2009) study, three possible rupture areas were considered for each modelled earthquake using the empirical scaling relations of Wells and Coppersmith (1994), which is not a subduction zone-specific relationship.
Following the study of Burbidge et al. (2009), several articles have been conducted on the PTHA in the MSZ. Using a limited mixed earthquake catalogue of the Makran region including shallow crustal, interface and intraplate events to constrain the truncated Gutenberg-Richter (TGR) recurrence relation, Heidarzadeh and Kijko (2011) defined three scenarios with identical fixed magnitude (Mw 8.1) in the eastern, central and western Makran for a probabilistic tsunami hazard assessment along the coasts of Iran, Oman and Pakistan. They considered fixed magnitude, location, and depth and also uniform slip distribution for each scenario.  Løvholt et al. (2012) conducted a tsunami hazard assessment on the global scale including the MSZ. Hoechner et al. (2016) used the TGR recurrence model and Makran instrumental era catalogue in computing the exceedance rate of tsunami height for Makran coastlines. The slip rate of the MSZ was ignored in their analysis, and they used the general seismicity of the whole Makran region without any distinction between shallow crustal, interface and intraplate earthquakes. El-Hussain et al. (2016 constructed some source models with different magnitudes along the MSZ and used the logic tree procedure and probabilistic method to estimate the tsunami hazard along the Oman shorelines. They only employed the TGR recurrence model for magnitude-frequency and assumed uniform slip on earthquake sources. Also, while not considering the coupling coefficient, probabilistic results were presented for the coasts of Oman. In another global-scale study, Davies et al. (2017) performed a probabilistic tsunami hazard assessment. They assigned three slip rates of 5.4, 9.0 and 13.0 mm/year for the Makran fault. Assuming a convergence rate of 20 mm/year, the resulting coupling factors are about 0.3, 0.5 and 0.7. They also assumed three maximum magnitudes of 8.1, 8.8 and 9.5. They did not take into account the segmentation of the MSZ into the western and eastern segments.
Another probabilistic tsunami hazard assessment was done by Rashidi et al. (2020b) for a range of heterogeneous slip distributions of the MSZ. Tsunami hazard curves and probability maps for the region were their outcome. They used the TGR recurrence model and a fixed non-planar geometry for the MSZ. The effect of seismic coupling was not a part of their computations.
In general, probabilistic tsunami hazard assessment is a multidisciplinary problem which contains tsunami numerical modelling and probabilistic calculations for the earthquake source and tsunami wave amplitude. This would require different experts in various fields of seismology, seismotectonics, geodesy, probability and statistics, hazard assessment, and numerical hydrodynamic modelling. So far, a comprehensive study, using the knowledge of the above various fields and combining them together in a logical and scientific format, has not been done at a regional-scale for the western Makran coastline. Present study tries to resolve and address the weaknesses and limitations of the previous studies as much as possible.
In this study, we consider the uncertainty of maximum seismic magnitude, earthquake occurrence models, continuity of seismic zone, seismic coupling coefficient, rupture depth, presence or absence of splay faults, fault locations, and fault slip distribution in PTHA for the MSZ. To clarify the importance of considering these uncertainties, it is worth noting that in modelling the tsunami observed after the only instrumental era large interface event in this region, i.e. the 1945 Makran earthquake with a magnitude of 8.1, there are significant differences between the hypotheses of various studies (e.g. Heidarzadeh et al. 2008Heidarzadeh et al. , 2009Jaiswal et al. 2009;Neetu et al. 2011;Heidarzadeh and Satake 2015).Considering homogeneous slip model, and presence or no presence of splay faulting (reaching the rupture to the seabed), are among the main reasons for the failure to model the observed tsunami wave height of this earthquake. This would place another emphasis on the key role of nonhomogeneous slip distribution to have a more reliable estimate of expected tsunami wave amplitudes, a factor that was ignored in most of the past PTHA studies (Table 1). In addition to non-uniform slip, the present study incorporates multiple assumptions including different top depths of the rupture, fault rupture-scaling model, the effect of coupling in seismic activity of the region, and a more precise bathymetry data (see Sects. 2). Earth subsidence effects are also evaluated in a probabilistic manner in this study for the first time.
It should be mentioned that our study only includes earthquake scenarios from the MSZ. As shown by the simulation results in Okal and Synolakis (2008), far-field scenarios have only mild effects on the Makran coasts and are thus excluded from our study. Moreover, the present study is only devoted to tectonic tsunamis, and thus, other types of tsunamis such as those generated by landslides are not considered. It is noted that the Makran region is also susceptible to landslide tsunamis as previously reported by Satake (2014, 2017).

Uncertainties in western Makran tsunami modelling
There are both aleatory and epistemic uncertainties in a PTHA study, though sometimes the margin between these two categories is vague and depends on personal opinion.
Aleatory uncertainties in fact indicate the natural randomness inherent in the nature of phenomena, whereas epistemic uncertainties are due to a lack of information and human knowledge about the natural processes. The logic tree is an efficient tool for considering epistemic uncertainties in probabilistic hazard and risk analysis studies, first introduced to the field by Kulkarni et al. (1984). One of the most important parameters in determining the behaviour and output of the logic tree is the weighting done by experts into its various branches, which represent different hypotheses in modelling. The branch weights in a logic tree represent the degree of belief of experts in each element of models.
In this section, by expressing the uncertainties in tsunami hazard analysis in the MSZ; the considered logic tree for this study Fig. 2 is introduced. It should be noted that, the dip angle of the subduction interface of the MSZ could be another source of uncertainty and was a matter of debate for a long time (see, e.g. Byrne et al. 1992, Maggi et al. 2000, Kopp et al. 2000Penney et al. 2017;Heidarzadeh et al. 2009). Penny et al. (2017 implemented well-located earthquakes in the region, and receiver functions from local seismometers to constrain the range of possible interface dips. They reported maximum average dips of ∼11• in the western, ∼9• in the central and ∼8• in the eastern Makran region, respectively. In a recent study, using the wide-angle reflection technique along three crustal-scale, trench-perpendicular, deep seismic sounding profiles, Haberland et al. (2021) has delineated the position and the dip of the subducting plate (oceanic Moho). According to Haberland et al. (2021), the crustal structure and the very gentle dip (∼8° ± 2°) of the subducting plate suggest a very wide contact zone, potentially allowing very wide asperities. In the light of aforementioned studies, here a fixed value is used for the dip of the subduction interface which is important for assessing tsunami wave height generated by different seismic scenarios.

Maximum seismic magnitude of the MSZ
For years, seismologists believed that some subduction zones would never be able to produce earthquakes of magnitude Mw 9.0 or greater (Ruff and Kanamori 1980). However, after the 2004 Indonesia Mw 9.2 and the 2011 Japan Mw 9.0 earthquakes, the general view of seismologists is that all or at least most of the subduction zones are capable to produce an earthquake of magnitude Mw 9.0, in a global analogy context. Generally, there are two main approaches for determination of the maximum magnitude (Mmax): those that 1 3 are based on an earthquake catalogue (maximum observed magnitude plus an increment) and those based on empirical scaling relations through a seismic source zone model. The second method, in the simplest way, uses fault length, while a more accurate alternative is to use the rupture area (which is proportional to the accumulated moment, Hanks and Kanamori 1979). Due to lack of historical events, except the 1945 Mw 8.2 earthquake in the eastern Makran, here we emphasize more on the second approach, which implies adopting the concept of global analogy of Mmax for all subduction zones (Meletti et al. 2010). As the seismic moment released in an earthquake is a function of rupture area, so assuming an average shear modulus of μ = 30 GPa for subduction zones (Hanks and Kanamori 1979), average dip angle of 8° and length of 900 km for the entire MSZ (or even the eastern or western segment of the MSZ), occurrence of an earthquake magnitude equal to or greater than Mw 9.0 is possible, as also previously reported by Heidarzadeh et al. (2009) Moreover, thermal modelling performed by Smith et al. (2012) to determine areas with temperatures above 150˚C in sediments showed that the occurrence of earthquakes with magnitude up to Mw 9.2 in the MSZ is possible. According to the Smith et al. (2012) modelling, the Makran crust reaches a temperature of 350˚C at a depth of 60 km and a temperature of 450˚C at a depth of 75 km. Therefore, with a length of about 800 km and a rupture width of about 200 ~ 350 km, a magnitude Mw 9.0 to 9.2 is proposed as the maximum Fig. 2 The logic tree used in this study magnitude of the MSZ in the worst-case scenario. In another study, by modelling deformation based on the GPS data to determine the downdip locking width, Frohling and Szeliga (2016a, b) proposed a maximum magnitude of Mw 8.8 for the entire length of the MSZ. Also, considering the length of more than 900 km for the MSZ and using empirical scaling relationships specific to subduction zones, such as that proposed by Skarlatoudis et al. (2016), occurrence of earthquakes with a magnitude of about Mw 9.1 ~ 9.2 is probable in the MSZ.
On the other hand, study of the location of megathrust subduction earthquakes shows that earthquakes with a magnitude of Mw 9.0 and more occurred in places such as the MSZ that have previously experienced quiescent periods with few numbers of moderate earthquakes (Heuret et al. 2011).
In a recent two-dimensional numerical study of megathrust earthquakes in subductions (magnitude greater than Mw 8.5), Muldashev and Sobolev (2020) showed that the dip angle of the subduction plate has the greatest effect on the maximum expected magnitude of a subduction zone. Hence, the lower the dip angles of the subduction zone, the greater potential for producing larger megathrust earthquakes. Therefore, taking into account low dip angle (~ 8˚) of the subduction plate in the MSZ, there is potential for occurrence of larger megathrust earthquakes.
Many earthquakes have been reported compared to smaller ones. Despite the extensive data, attempts have been made to obtain empirical scaling relationships between fault rupture dimensions and seismic magnitude based on past earthquakes data (Papazachos et al. 2004;Blaser et al. 2010;Strasser et al. 2010;Murotani et al. 2013;Skarlatoudis et al. 2016;Allen and Hayes 2017;Thingbaijam et al. 2017). One of the important points is whether the scaling relations, which are mainly constrained and developed by using smaller earthquake data, can also be attributed to megathrust earthquakes in subduction zones. Actually, there is still low data to determine the rupture behaviour of giant earthquakes (magnitude greater than 8.5). This limitation is especially evident in determining the width of the rupture plate (Skarlatoudis et al. 2016). Some researchers suggest rupture width saturation, meaning that the fault width does not change by increasing the earthquake magnitude after a certain value (Tajima et al. 2013;Skarlatoudis et al. 2016). However, the lack of data in this range of magnitude makes it difficult to draw final conclusions.
According to Strasser et al. (2010), it is important to consider more than one empirical scaling relationship to determine maximum magnitude from an epistemic uncertainty perspective. Therefore, in this study, the possible maximum earthquake magnitude has been calculated with three scaling relationships of Skarlatoudis et al. (2016), Thingbaijam et al. (2017) and Allen and Hayes (2017), based on different assumptions regarding the fault rupture length, width and area.
These three selected relations have considered more earthquakes than the other existing scaling relations, including very large recent earthquakes (for example, the 2011 Tohoku earthquake) and therefore are more reliable in this respect.
In this study, the selected scaling relationships are used with equal weight in determining the maximum magnitude for the MSZ based on the width, length and area of the rupture plate as follows: • Empirical scaling relationships based on magnitude-rupture area (considering two rupture area scenarios). • Empirical scaling relationships based on magnitude-rupture width (considering two rupture width scenarios). • Empirical scaling relationships based on magnitude-rupture length.
Therefore, if the rupture width (width of the locked area) is assumed between 80 km (minimum depth of 14 km and maximum depth of about 25 km) and 185 km (minimum depth of 14 km and maximum depth of about 40-45 km), based on the above-mentioned empirical scaling relationships, the maximum possible magnitude for the two pessimistic and optimistic scenarios in each of the western, eastern and entire Makran logic tree branches will be (8.7, 8.5, 9.0), (8.3, 8.3, 8.5), respectively, for the central magnitude in the pure characteristic model. The maximum moment magnitude in the GR models for each of the western, eastern, and entire Makran logic tree branches are selected equal to (8.7, 8.7, 9.1), (8.3, 8.3, 8.5), respectively, for pessimistic and optimistic scenarios. In this regard, the maximum moment magnitude in the GR models and the central magnitude in the pure characteristic model for each of the western, eastern, and entire Makran logic tree branches are selected equal to these values.

Earthquake occurrence model
In a PTHA study, an empirical model that gives the number of earthquakes with magnitude m that occur within a year is needed to explain the magnitude-frequency distribution of earthquakes. In earthquake hazard studies, either the GR or the characteristic earthquake models are commonly used. As pointed out by Davies et al. (2017), two modified versions of the GR models, that incorporate a maximum magnitude bounds, are often used to characterize earthquake recurrences: the model introduced by Kagan (2002) and the TGR model as follows: The above model is also known as the exponential model in the PTHA studies and has a zero recurrence rate at Mmax. The pure characteristic model, also known as the "Maximum magnitude" model, generates earthquakes only within a relatively narrow magnitude range (Stewart et al. 2002). As earthquakes of small-to-moderate magnitude that follow the GR distribution are unlikely to cause a tsunami (Thio et al. (2007)) in addition to the TGR model, we also use the pure characteristic model (Youngs and Coppersmith 1985).
The pure characteristic model is more consistent with the observed seismicity of the Makran region (lack of moderate to large events) and also does not require a b-value (see Sect. 2-3). This is important from the viewpoint that some recurrence intervals resulted from a GR model, such as one earthquake with M ≥ 7.5 every 30 years, which are reported by El-Hussain et al. (2016), cannot be validated at least through the last 200 years (Ambraseys and Melville 1982).
A pure characteristic recurrence model with a width parameter of 0.3-0.5 magnitude units (depending on the past seismicity) was used in a probabilistic tsunami hazard study in Japan, representing characteristic earthquakes with a uniform distribution (Annaka et al. 2007).
Because only earthquakes with magnitudes greater than 7.5 can generate tsunamis, Thio et al. (2007) in another study on the probabilistic hazard of tsunamis in eastern Asia rely solely on the pure characteristic models. A similar study is linked to what is known as quiescent seismicity in the Cascadia region of eastern Canada and USA, where the GR model was not used; however, the recurrence model for the interface earthquakes was recognized as a pure characteristic model with a return period of 500 years for a M > 9 event and a return period of 100 years for a M > 8 event (Frankel et al. 2002). Atkinson and Goda (2011) used an alternative version of pure characteristic recurrence model, assuming a normal distribution with a mean of 8.5 and a standard deviation of 0.5 to calculate the recurrence rates for various moment magnitudes of the Cascadia subduction zone.
To constrain the recurrence models using fault slip rates, the following relationship between the annual rate of moment-release and earthquake recurrence rates is suggested by Anderson and Luco (1983): where ̇n(m) is the density function for earthquake occurrence rate, Ṁ 0 denotes the annual rate of moment-release, and M 0 (m) denotes the magnitude-seismic moment relationship proposed by Hanks and Kanamori (1979): It is worth noting that both the GR and the characteristic models were favoured in different studies when deciding on a recurrence model (Schwartz 1999;Pisarenko and Sornette 2004).

b-value in the Gutenberg-Richter model
Following the discussion above, two recurrence models, namely the TGR and the pure characteristic model (also referred to as maximum magnitude recurrence model), are used in the current study. Therefore, to completely explain the seismicity in the Makran subduction region, we must obtain the b-value.
According to Benito et al. (2012), in any subduction zone, three seismicity regimes, crustal, interface slab, and inslab, can cause earthquakes. In the Makran region, the lack of interface slab earthquakes is a significant limitation to use the GR recurrence relation. Unfortunately, due to a lack of data, seismicity parameters in the subduction zone are normally calculated using a combination of interface, inslab and crustal earthquakes. For example, Deif and El-Hussain (2012) and Hoechner et al (2016) combined earthquakes from all three seismic regimes and estimated the b-value to be 0.77 and 0.82, respectively. Also, a b-value of 0.6 was calculated by Heidarzadeh et al. (2008) and Heidarzadeh and Kijko (2011) for the Makran subduction zone, although without distinguishing inter and intraslab earthquakes.
Another issue caused by the lack of earthquakes, as emphasized by Pisarenko and Sornette (2004), is weak constraint at the end tail of the GR distribution, where major tsunamigenic events occur. It is worth noting that the aforementioned issues also exist in other subduction zones, such as the Cascadia subduction zone, where seismic activity is currently low.
Since the lack of seismic activity prevents the calculation of the b-value in the Cascadia area, in a comprehensive study by BC Hydro (2012), its value was adopted from other studies around the world where the abundance of seismicity permits the calculation of the b-value. The absence of interface seismic data and the presence of a thick sediment layer According to the GEM Faulted Earth Subduction Interface Characterization Project: Version 2.0 (Berryman et al. 2015), who processed seismicity data from subduction zones around the world and took into account the low rate of M > 5.5 earthquakes in Hikurangi, the Caribbean, Cascadia, and Makran, they suggested that the b-value for the interface slab in these regions is likely lower than the global average.
By reviewing the b-values reported worldwide, Wang et al. (2016) allocated 0.7 to the north Taiwan subduction zone, Power (2013) used 0.75 for the Japan Trench, and GEM group (Berryman et al. 2015) used a range of 0.7-1.2 for subduction zones around the world. Note that, in the GEM report, the b-value in Makran was assigned according to Bayrak et al. (2002) who devised the b-value equal to 0.78. However, the combination of Zagros (a shallow crustal environment) and Makran earthquakes, which was disregarded in the GEM report, is a serious criticism of Bayrak et al. (2002).
We calculated the b-value by reviewing previous global studies and gathering seismic activity in the western and eastern MSZs. We obtained a b-value of 0.75 without distinguishing between shallow crustal and interface slab regimes, but excluding inslab events (Rashidi et al. 2020a).

Continuity of seismic zone
The MSZ is an area with an approximate length of 1000 km which has been extended through the coast of Iran and Pakistan from Strait of Hormuz to Karachi in Pakistan. According to (Heuret et al. 2011), when the length of subduction is below ~ 500 km, the large earthquakes are unlikely to happen, while longer subduction zones may potentially be able to extend the rupture along their whole length. However, the amount of rupture observed in past longer subduction earthquakes around the world has also been limited to a part of the total length of the subduction zone (500 to 2000 km). The reason for this discontinuity in rupture propagation, and consequently the limited rupture length, is that the complexity of the subducting plate (such as different sediment thicknesses, fracture zones, Seamounts, the presence of horst and Graben) and the strength of the upper plate cause local changes in friction along the length.
Although there are strong evidences of the segmentation of the MSZ into two parts, eastern and western, with different geodynamic behaviours (Byrne et al. 1992), the possibility of a very large earthquake due to the cascading rupture of both segments could not be ignored, assuming the rupture spreads from the west to the east or vice versa. Therefore, it is necessary to consider the MSZ as an integrated seismic region, although not very likely, in order to consider all possible scenarios for thorough tsunami hazard estimation.
In previous research works conducted, such as the tsunami hazard analysis for the coasts of Oman ( El-Hussain et al. 2016), an integrated rupture in Makran with a magnitude of more than 9.0 with a branch weight of 0.15 has been considered. Also, in the seismic hazard analysis report of the cities of Abu Dhabi, Dubai and Ras-al-Khaimah in the United Arab Emirates (Aldama-Bustos et al. 2009), a weight of 0.05 is attributed to the scenario corresponding to the cascading rupture of the MSZ.
In present study, taking into account the above discussion, a logic tree branch with a low weight (0.15) is considered to account for the possibility of a cascading rupture along two MSZ segments (the unsegmented branch in the logic tree of Fig. 2).

Seismic coupling coefficient
The seismic coupling coefficient s is the ratio of seismic moment released by earthquakes to seismic moment stored between plates as a result of relative tectonic movement within a subduction zone.
The seismic moment rate is generally calculated by dividing the total seismic moment obtained from recorded earthquakes in an instrumental catalogue by the catalogue's time period. It is worth noting that the estimation of seismic moment rate is inherently uncertain due to the long return periods of major earthquakes that might not be reported in instrumental catalogues. Actually, a century of seismic activity in a subduction zone does not adequately represent the long-term seismicity, and the approximation based on it, results in lower coupling coefficients. Taking into account the scarcity of seismicity in the Makran region, seismic coupling of the MSZ is another open question.
The measurement of the inter-seismic strain accumulation rate, which gives the interseismic moment rate ( X g ), has been easily accomplished due to the rise in the number of GPS stations (Scholz and Campos (2012)). The shortcomings of the traditional catalogue-based method for calculating the seismic coupling rate, which is primarily due to the lack of large earthquakes recorded in seismic catalogues, have been resolved by the novel geodetic approach. Note that, although the coupling ratio calculated by the classical approach yields the released moment, the geodetic coupling rate gives the total energy stored within the tectonic plates. Both methods are equal in the long run if the coupling rate remains constant over time. The inconsistency between the two concepts discussed above highlights earthquakes that have been ignored in seismic catalogues and should be explored further using historical resources or paleo-seismic methods (Scholz and Campos (2012)). One source of uncertainty in the estimation of coupling rate is the error in the approximation of coupled surface area A c , which is the region of the contacting plates where seismic slip occurs (seismogenic area). The length of this area L c is inferred from the length of the subducting segment. However, the end pieces from both sides are often ignored. A more important source of uncertainty is the estimation of the coupled area's width W c . W c is determined between the accretionary wedge and the brittle-plastic transition zone. According to Scholz and Campos (2012), the brittle-plastic transition normally occurs where the temperature reaches 350 °C. Another transition zone exists below this temperature, which is where slow-slip events usually occur. Further below, the plates can freely slip without any coupling.
There are three ways to calculate the width of the coupled field, W c : first, by observing the location of background seismic activity, second, based on the slipped areas of large previous earthquakes, and third, by inverting the GPS data. The first two methods have a 30 per cent uncertainty (Scholz and Campos 2012), while the third method has a lower uncertainty and is more effective in determining the deep seismic edge. Because the top seismic edge is usually located far from the coast and the GPS stations are located along the coast, the data provided by the GPS stations do not provide a good constraint for the top edge depth. Scholz and Campos (2012) used both GPS data and seismological observations to calculate the amount of seismic coupling in many subduction zones around the world. According to their study, the coupling rate is between 0 and 1, with an average value of 0.6.
In the other study, Heuret et al. (2011) reported the coupling coefficient based only on the previous century seismicity, and the slip rate values across the subducting area. It is worth noting that Heuret et al. (2011) used a shear modulus of G = 50 GPa instead of the more common values of G = 30 GPa (Hanks and Kanamori 1979;Davies et al. 2017) or G = 40 GPa (Scholz and Campos 2012), which resulted in a 20% reduction in the coupling coefficient values compared to what was reported by Scholz and Campos (2012). The seismic coupling by Heuret et al. (2011) was found to be uncorrelated (R = 0.17) with the subduction velocity, which is consistent with the observation of Pacheco et al. (1993).
In a comprehensive study by GEM group (Berryman et al. 2015), the coupling ratio in the MSZ was assumed to be between 0.3 and 0.7. Note that, in (Berryman et al. 2015), the coupling ratio was based on expert opinion rather than modelling. As emphasized by Byrne et al. (1992), because the creeping regime exists at every subduction zone, we assumed two low and high creeping branches for both eastern and western Makran. Furthermore, the fact that some faults have a cycle of seismic activity and silence was taken into consideration when assessing branch weights. Davies et al. (2017) considered subduction rate as 19 mm/yr and used three different seismic slip rates of 5.4, 9.0, and 13.0 mm/yr, corresponding to coupling rates of 0.3, 0.5, and 0.7, respectively. Furthermore, maximum magnitudes of 8.1, 8.8, and 9.5 were assumed, accordingly. Ghadimi et al. (2021Ghadimi et al. ( , 2022 computed the long-term crustal flow and forecast of seismicity for the MSZ, using a kinematic finite element technique. They investigated the possibility of creep in western Makran, by comparing two different models with the coupling rates of 0 (means free slipping) and 1 (corresponds to a totally locked fault) for the western Makran. The results of the temporarily locked model indicated about 1 ~ 3 mm/yr higher rate of shortening than the other models with the steady creeping subduction. Also, the predicted inter-seismic velocities and seismicity from the creeping model are more accurate for the western Makran.
In this article, based on the results of Ghadimi et al. (2021, submitted) and the review presented above, we used a coupling ratio of 0.3 for western, 0.6 for eastern, and two values of 0.3 and 0.45 for an unsegmented model. Additional conservative branches are also included in the logic tree for the eastern, western, and unsegmented model, with a pessimistic coupling coefficient of 1.0, indicating a fully coupled subduction (Fig. 2).

Updip rupture depth
The shallow aseismic-seismic transition depth ( d s ) estimated by Pacheco et al. (1993) as d s = 10 km for all subduction zones was adopted by the Slab1.0 Model (Hayes et al. 2012). Slab 2.0. (Hayes et al. 2018) takes a different approach based on the seismicity depth of subduction zones, with a minimum d s = 9 km (Tonga subduction zone), a maximum ds = 16 km (Kamchatka subduction zone), and an average ds = 12 km. In contrast to other active subduction zones around the globe (Smith et al. 2012), the MSZ has a thick sediment layer (3-7 km) (McAdoo et al. 2004), so an active depth deeper than the global average seems reasonable. For the MSZ, considering the limited seismic data obtained in the subduction interface, Slab 2.0 (Hayes et al. 2018) reports shallow and deep seismogenic depths of 13 km and 40 km, respectively. Taking into account the above discussion and considering depth distribution of local seismicity in the region (Rashidi et al. 2020a) (Fig. 3), the updip rupture depth in Makran has been set at 14 km as an optimistic value.
Prior to the Indian Ocean Tsunami of 2004 and the Tohoko Tsunami of 2011, the possibility of an earthquake rupture extending up to the seabed was thought to be unlikely. However, clear evidence suggests that the rupture was extended up, nearly to the subduction Fig. 3 Top: Focal mechanism and depth distribution of Makran seismicity. Bottom: Section of Makran seismicity, based on catalogues with more reliable depths including, normal faults, reverse faults, strike-slip faults, and without mechanism. A line is drawn with a dip angle of 8 degrees from latitude 24 degrees north, which indicates the boundary of the oceanic lithosphere of the sinking Arabian plate below Eurasia. Depths of 14 and 40 km are marked by two continuous horizontal lines. The ratio of vertical scale to horizontal scale is 1 to 3. The triangles show the volcanoes of Bazhman, Taftan, and Sultan in both figures trench in both cases, and was shallower than what previously thought. Although the rupture is unlikely to nucleate in shallow depths, it is now widely accepted that deeper events could extend seaward (Smith et al. 2012). Accordingly, the probabilistic study of tsunami hazards must account for seabed rupture. Therefore, in addition to the optimistic value of 14-km, we consider 3-km as an alternate updip rupture depth for buried rupture scenarios. In fact, the 14-km updip rupture depth assumes that the rupture would not spread into unconsolidated shallow sediments; whereas the 3-km ones takes this into account.
We also consider trench faulting as another possible option, suggesting the rupture could extend up to the seabed. Splay faulting is also considered in this study according to Sect. 2_7 (for different models of faulting means buried, splay and trench-breaching see Gao et al. 2018). Thermal and geodetic constraints can help determine the deep extent of the seismogenic zone, but expert debates are generally determining the shallow extent (e.g. Stirling et al. 2012 for Hikurangi subduction interface).

Splay faulting
Splay faults are secondary faults rising from the plate interface and extend up to the seafloor. They are usually detected in seismic reflection profiles and exist in many of the subduction zones, including Makran (Grando and McClay 2007;Smith et al. 2012;Mokhtari 2015;Heidarzadeh and Satake 2017). Splay faults can slip during large subducting events or independently produce earthquakes. They are usually shallow (or even reach the surface) and steeper than the megathrust and can locally and significantly displace the seafloor and therefore the initial water surface. Subsequently, they can have a substantial contribution in tsunami hazard at local scale. The 1946 Nankai tsunami is a good exemplar of splay fault tsunamis (Cummins and Kaneda 2000). The effect of splay faulting in enhancing the 1960 Chile, 1964 Alaska and 2004 Indian Ocean tsunamis has been also reported (e.g. Plafker 1972;Sibuet et al. 2007). Splay faulting was brought up as a secondary tsunami source that increased the run-up caused by the 1945 Makran tsunami (Heidarzadeh et al. 2009). The dip of a splay fault can be landward or seaward in which they both can enhance tsunami wave amplitudes. However, a seaward dipping fault can generate a stronger tsunami as the shallower part of the fault is towards the sea. Priest et al. (2010) took into account the effect of splay faulting in evaluating tsunami hazard of local earthquakes in the Cascadia subduction zone. To do so, they developed a logic tree and assigned different weights for their scenarios, including splay faulting. In the case of splay fault scenario, the slip is partitioned between the megathrust and splay fault. Priest et al. (2010) showed that the splay faulting amplifies the tsunami wave height between 6 and 31% and the run-up between 2 and 20%. The results of tsunami numerical simulation for different models of faulting (buried, splay and trench-breaching) in Northernmost Cascadia (Gao et al. 2018) show that the splay faulting model generates higher tsunami waves (50 to 100%) than the commonly used model of buried rupture. Heidarzadeh et al. (2009) studied the effect of splay faulting on tsunami amplitudes in Makran. Their outcome for a combined tsunami source of the eastern Makran segment and a splay fault indicates the increase in seabed displacement by 1.5 times and tsunami amplitude by two times at some locations.
In this study, we incorporate the effect of splay faulting into the tsunami hazard along the coastline of Iran. We consider a landward splay fault dip of average 35° (ranging from 20˚ to 55˚ for different ruptures, which is reasonable according to Grando and McClay 2007;Smith et al. 2012) that is joined to the megathrust (main fault) at a depth of about 1 3 14 km (Fig. 4). We partition the seismic moment of the rupture between the main fault and the splay fault in which the splay faulting portion ranges from 10% for the largest scenarios to 32% in the case of smallest scenarios.

Fault rupture locations and slip distributions
Since the occurrence of an earthquake with different magnitudes, in any parts of the MSZ is possible, different fault locations along the subduction interface should be considered to cover the entire seismic zone. In this way, more than 155 fault locations which completely cover the MSZ are assumed in this study (Fig. 5).
Previous studies (Goda et al. 2014;Ruiz et al. 2015;Sun et al. 2018;Crempien et al. 2020;Momeni et al. 2020) have shown the amplifying effect of heterogeneous slip distribution on the generated tsunami wave height. In order to model the heterogeneous slip on the fault, the finite fault technique (Gallovič and Brokešová 2004) has been used The lengths and widths of the rupture plane for each earthquake scenario and the mean and maximum slip values are calculated using appropriate empirical scaling relationships (Allen and Hayes 2017).
To consider aleatory uncertainty of slip distribution, for each fault location, 12 heterogeneous slip scenarios are generated using k-2 kinematic source modelling (Gallovič and Brokešová 2004), as according to the Tsunami Working Group in (2006), considering at least 12 slip scenarios is an appropriate and reliable choice. The heterogeneous slip scenarios are generated considering limitations of average and maximum slipping on fault which are calculated using Allen and Hayes (2017) scaling relationships. Also, taking into account the previous inversion-based earthquake rupture models of interface events (Mai and Thingbaijam 2014), in generating heterogeneous slip scenarios, placement of slip asperities in the upper one third of the fault in more than 3 of 12 scenarios is not allowed. In other words, in most of scenarios, slip asperities are located at the deeper parts of the fault. Figure 6 shows six samples of generated slip distribution scenarios in this study corresponding to Mw 8.0, Mw 8.4, and Mw 9.0.

Partial summery
In this section, we expressed the uncertainties in tsunami hazard analysis in the MSZ and introduced the logic tree considered in this study (Fig. 2). Due to the uncertainties of the seismic magnitude, fault location, slip distribution, seismicity depth, and splay fault, we have considered a number of scenarios according to Table 2 that perform numerical simulations (Sect. 3).

Tsunami numerical modelling
The tsunami phenomenon consists of three main phases: generation, propagation, and inundation. Underwater earthquakes are the main cause of about 80% of tsunamis (Løvholt 2017). Therefore, we consider only these kind of earthquakes as tsunami sources.
If the total rupture time of an earthquake is much shorter than the period of the tsunami wave, the seabed displacement is considered as an instantaneous motion. This method is conventionally called passive seismic generation (Dutykh et al. 2006;Kervella et al. 2007). By using this method, static (rather than dynamic) displacement is considered for the seabed due to earthquake. It is also assumed that the initial amplitude of the tsunami wave is exactly the same as the amount of vertical displacement of the seabed, because the water of this vast area affected by fault sliding cannot be drained within the rupture time (Steketee 1958). Therefore, the initial conditions for the propagation of tsunami waves are obtained by transferring the final vertical static displacement of the seabed to the water surface. Okada's (1992) closed form solution is used to obtain vertical displacement in the seabed. According to this solution, in a semi-infinite homogeneous environment, if the fault is rectangular and the amount and direction of slip are constant everywhere on the fault, a closed form solution can be written to obtain displacement anywhere in the space. For the faults with non-constant slip, the fault is divided into a number of sub-faults, that can be treated separately and the final displacement of the seabed is equal to the sum of all sub-fault displacements based on the superposition principle.
The source rupture modelling parameters used in this study are listed in Table 3. After tsunami waves are generated, they propagate through the water. The most common equations to describe the propagation of tsunami waves are the shallow water equations, which are derived from mass and momentum continuity. These equations assume that the horizontal scale (wavelength) is much larger than the vertical scale (water depth). The second assumption for these equations is that the viscous effect is ignored. As a result of these two assumptions, the horizontal velocity of the waves will be much higher than the vertical velocity of the waves, the horizontal velocity of the wave will be constant at depth and the pressure will be hydrostatic. This study adopts COMCOT (COrnell Multi-grid COupled Tsunami Model; Liu et al. 1998) as the kernel to perform tsunami modelling. COMCOT supports both linear/nonlinear shallow water equations for tsunami simulations with the two-grid nesting scheme for flexible computational grid sizes and time steps crossing multiple layers and moving boundary scheme for calculating coastal inundation areas and flood depths. On the one hand, in deep-water regions where the nonlinear effect is ignored, the linear shallow water equations are used for simulating a proposed tsunami from its original source to nearshore regions. On the other hand, in shallow water regions where tsunami waves are amplified by the shoaling effect and bathymetry effect (see, e.g. Imamura et al. 1988;Kajiura and Shuto 1990), the nonlinear equations are adopted for calculating nearshore tsunamis. It is noted here that both linear and nonlinear shallow water equations are conducted in the spherical coordinate.
As the wavelength of tsunami waves changes due to the water depth, the grid sizes and time steps shall be modified accordingly. Thus, we consider the three-layer nested grid configuration to simulate tsunami waves propagating from offshore to nearshore (Fig. 7).
The grid sizes of the first, second, and third layers are 30, 15, and 5 arc-seconds, respectively. The bathymetry data for the first and second layers come from the General Bathymetric Charts of the Oceans (GEBCO) (2020) digital grid atlas (IOC 2003;Weatherall et al. 2015), which has the spatial resolution of 15 arc-seconds. The data used in the third layer come from local data in the Chabahar area with a high spatial resolution of 3 arc-seconds. It should be mentioned that we use the Manning's coefficient equal to 0.025 as a standard value for coastal water nearshore bottom friction in layer 2 and 3. The computational settings of the nested grid configuration can be found in Table 4. 1 3 The calculation of the water amplitudes at the shoreline using low-resolution bathymetry data may underestimate the tsunami impact (Kamigaichi 2009(Kamigaichi , 2015. Since in most of the coastal areas of western Makran (except the Chabahar Bay) bathymetry data are only available with a maximum resolution of 15 s, the Green's law is applied in this study, ensuring the calculated shoreline results would not be less than the actual expected values. Moreover, for a regional study including a large number of widespread locations distributed over hundreds of kilometres, refined numerical inundation simulations are too time-consuming to perform for all scenarios and locations (Davies et al. (2017), andLøvholt et al. (2012)). Løvholt et al. (2012) in a comprehensive study implemented a procedure to transfer the surface elevations to the maximum shoreline water levels by determining a set of amplification factors based on the parameters for the incoming wave and some idealized bathymetric profiles. Their results show a good estimation of coastline water height using the Green law's method comparing to the numerical methods with different nearshore bathymetric profiles. Moreover, whereas considering a no-flux boundary condition at the shoreline for amplification factor method leads to a double value of surface elevation due to reflection, Green's law results are good and reliable.
According to the Løvholt et al. (2012) study and taking into account the typical bathymetry of the Makran coastal region, it could be argued that in most locations, using the Green's law for estimation of coastal line water height would lead to reasonable results. In addition, we note that the Green's law method has also been adopted to   Hoechner et al. 2016;Davies et al. 2017).
According to the Green's law (conservation of energy), the tsunami wave height farther from the shore is used to approximate the tsunami wave height on the shore according to Eq. (2). h and d denote the tsunami wave height and see depth, and indices '0' and '1', respectively, represent the value on the shoreline and offshore point.
In this study, the height of water on the shoreline is calculated by converted values from depth of 50 and depth of 10 m by using Green's coefficients. The final result for coastal line is presented using these two values with reasonable weights.
A number of limited scenarios were considered to select the time required to simulate. The modelling results show that the maximum amplitude of seismic waves for all stations occurs before the first 6 h and the simulation time of 6 h is a good choice. Fig. 8 The maximum water height calculated on the coastline, which is calculated by converted values from a depth of 10 m for all numerical simulations (1776 scenarios)

3
The results of numerical simulations (1776 scenarios) for coastline, which is calculated by converted values from depth of 10 m by using Green's coefficient, are shown in Fig. 8.
After completing the numerical simulations for all scenarios, the PTHA calculation is done according to Sect. 4.

Probabilistic tsunami hazard analysis (PTHA)
Considering the uncertainties of the seismic magnitude, fault location, slip distribution, seismicity depth, and splay fault, we have considered 1776 scenarios and have performed numerical simulations for them according to Sect. 3. Based on the results of the numerical simulations, for each magnitude m and each wave height z 0 , the probability of the wave height exceeding a given value of z 0 , P Z > z 0 |m can be obtained. This probability is obtained by integration over the aleatory uncertainties representing the natural uncertainty in the expected wave amplitudes given a specific magnitude of tsunamigenic earthquake. Hazard curves are obtained for different branches of logic tree representing epistemic uncertainty.
For each branch of the logic tree, the annual exceedance rate ν(z 0 ) is calculated based on the following equation.
where (m) is the annual frequency of earthquakes with magnitudes m ± Δm∕2. The magnitude bin size Δm is set as 0.2 magnitude units (see Table 3).
In the logic tree branches where the activities of the eastern and western zones of Makran are considered different, the random slip distributions, related seabed deformations, and resulting tsunamis for each zone are simulated separately during calculation of the exceedance rate. For these branches, the exceedance rate ν(z 0 ) is calculated based on the following: where N Segments is equal to 2 (eastern and western MSZs) and n (m) is the annual frequency of earthquakes on source n with magnitudes between m ± Δm∕2.
The calculation of P (Z > z 0 | m) can be done in two ways. The first method is to form a set of simulated wave heights at the desired point, which are obtained by considering all scenarios with in the magnitude bin m ± Δm∕2 . Then, a wave height like z 0 is selected and the number of wave heights that is greater than z 0 is counted in the set. The ratio of the number of scenarios with wave heights greater than z 0 to the total number of the scenarios provides the intended probability at the desired point, which means as the relative frequency based on the simulated samples.
In the second method, it is assumed that the desired probability follows the log-normal distribution (e.g. Kulkarni et al. 2016). Therefore, by calculating the mean and standard deviation of the set of simulated wave heights, the log-normal distribution is fitted to them. Figure 9 shows the comparison of the relative frequency based on the simulated samples distribution and log-normal distribution in calculating P (Z > z 0 │m), which is presented Fig. 9 Comparison of real distribution and log-normal for a point on the coast of Chabahar according to the simulation results for a point on the coast of Chabahar and shows that two approaches give similar results, at least for probabilities higher than 0.05. Then, based on the assumptions of the logic tree, we calculated the probabilistic tsunami hazard curves for all of its branches and formed the statistical population for different hypotheses of the exceedance rate. Based on the formed sample population, for each wave height, we calculated exceedance rate at the average, 16%, 50% and 84% percentiles of distribution. Figure 10 shows an example of these calculations for a point located on the coast of Chabahar. By reading the wave height in the values of 1/475 and 1/2475 from the output hazard curve, it is possible to calculate the value of the wave height in the 475-year and 2475-year return periods, at different levels of confidence. It is worth noting that a group of hazard curves relating to the eastern Makran scenarios are located below others in Fig. 10. The reason for the lower hazard is that Chabahar is in the western Makran, whereas all of the source zones in those scenarios are located in the eastern Makran, therefore, the sourcereceiver distance is comparatively longer for this group. Fig. 10 An example of calculating exceedance rate at the 16%, 50% and 84% levels of distribution for a point located on the coast of Chabahar, for the 475 and 2475 years return period. Note that, a group of hazard curves relating to the eastern Makran scenarios are stacked below the others 1 3

Earth subsidence
Here, the seismic source is a local subduction event (as is possible in the MSZ), in addition to the uplift of the sea floor, it can cause land subsidence in areas near the coast or on the coast (ASCE7 2016; Ruiz et al. 2015). To take into account this effect, the change in the ground elevation per point source is calculated using the Okada (1992) method and the superposition principle is employed for extended fault ruptures. The subsidence is added to the tsunami wave height on the shore and at the shoreline to retrieve the final results (ASCE7 2016).
In this article, in order to have more realistic assessment of the tsunami water height relative to the reference ground level, we perform probabilistic hazard calculations once for the height of water from the water free surface and once for the earth subsidence, and sum up the results of the two, conservatively. Actually, this is somehow conservative, as we know that two maximums do not occur simultaneously. In each modelling, if there was an uplift at any point along the shoreline instead of subsidence, the subsidence value was considered zero. The results of earth subsidence for western Makran coast, in the vicinity of Sistan va Baluchestan and Hormozgan provinces, have been calculated in two return periods of 475 and 2475 years with two levels of average confidence and eighty-four per cent (Figs. 11 and 12).
According to this calculation, the maximum subsidence of the coastline with average level of confidence is 0.15 m and 1 m for return periods of 475 and 2475 years, respectively.
Moreover, the maximum subsidence of the coastline with 84% level of confidence is 0.5 m and 1.23 m for return periods of 475 and 2475 years, respectively.

Results and discussion
In this study, the final results of tsunami wave height for western Makran coast have been calculated in the form of uniform hazard maps in two return periods of 475 and 2475 years with two confidence levels of mean and eighty-four per cent (Figs. 13 and 14).
The output presented in this study for each 15 km section of the Makran coast is equal to the maximum amount of probabilistic analysis observed in all points located in that section.
Since the return period of 2475 years is more compatible with the long recurrence times of large plate boundary earthquakes in the MSZ, it is recommended to use the tsunami wave heights of 2475 years (at the mean level) for structural design purposes. This is consistent with the "maximum considered tsunami" return period adopted by the ASCE7 (2016) for the US coastal regions. The more conservative level of 84% could be applied for critical structures.
Comparison of the output results of our simulations (before applying the probabilistic hazard analysis algorithm) with simulated scenarios with similar dip angle (8 to 15 degrees) in other places such as the scenarios of the 1707 Nankai earthquake in Japan (Furumura et al. 2011) and the coast of Sumatra, Indonesia (McCloskey et al. 2008) indicates the similarity of the range of obtained wave height. Moreover, the accuracy of the probabilistic hazard analysis process is examined and compared with the results for different coupling coefficients (lower pairing coefficient; lower wave height), different return period (lower return period; lower wave height) and different confidence levels (lower confidence level; lower wave height).
Finally, it is a good idea to compare the results with some previous regional or global studies, devoted to the PTHA in the MSZ. In comparisons, it should be kept in mind that the subsidence effects are only included in our results.
The probabilistic hazard maps reported by Burbidge et al. (2009) are divided into "high-hazard" and "low-hazard" maps based on different assumed maximum magnitudes. They noted that the actual hazard occurs between these two scenarios. The western segment of the MSZ is not included as a source zone in their low-hazard scenario, because of the fact that there is no certain knowledge of its potential to produce large tsunamigenic earthquakes. The highest tsunami amplitudes (in the water depth of 100 m) from these two low-and high-hazard scenarios along the coastline of Iran are 0.3 and 2.7 m. Thus, for a return period of 2000 years, applying the Green's law approximation (Kamigaichi 2015), the coastal tsunami amplitudes would be 1.0 and 8.5 m, respectively. In a more recent study on the PTHA of the Makran region, Hoechner et al. (2016) give tsunami wave heights of 1 and 3 m (for the return period of 500 years) and 4 and 11 m (for the return period of 2500 years), at the Jask and Chabahar cities that are comparable with our results.  Davies et al. (2017) conducted a PTHA on the global scale including the MSZ and showed that the probabilistic tsunami wave height along the Makran coastlines is generally about the same range, but it decreases towards the Strait of Hormuz. In their results, tsunami wave amplitude varies between 1.0 and 5.0 m and between 5.0 and 10.0 m for return periods of 500 and 2500 years, respectively, which is comparable with our results.

Conclusions
This study presents the first tsunami hazard zonation in the coastal area of western Makran using probabilistic tsunami hazard assessment, including all possible uncertainties and earth subsidence effects. The results are expressed in terms of uniform hazard maps for two return periods of 475 and 2475 years at two confidence levels of average and eighty-four per cent. As tsunamis generated by far-field tectonic sources cause very little effect on this part, we consider only tsunamis of earthquake origin generated in the MSZ.
The uncertainty of maximum seismic magnitude, scaling relations, earthquake occurrence model, continuity of seismic zone, seismic coupling coefficient, rupture depth, splay fault, fault locations, and fault slip distribution are included in PTHA calculation for this region. We use the COMCOT tsunami modelling code to simulate 1776 scenarios.
The result shows that Chabahar and Sirik are, respectively, exposed to the highest and lowest probable water heights due to the tsunami.
In this study, we do our best to have reliable assessment for PTHA in the west MSZ. However, it is important to consider its limitations and try to reduce them as possible in future studies. In this study, we consider underwater earthquakes as the only sources for tsunami, but it seems appropriate to study seabed landslides as another possible causes of tsunami and their combination with earthquakes in the MSZ as previously reported by Satake (2014, 2017). In addition, in our study, we ignore the dynamic interaction between the tides and tsunami waves. Therefore, to use the probabilistic tsunami water height presented in this study, it is necessary to add the value of corresponding maximum tidal height at any location, conservatively. More accurate results can be achieved by considering this dynamic interaction. Moreover, in this article, we perform probabilistic hazard calculations once for the height of water from the water free surface and once for the earth subsidence, and sum up the results of the two, conservatively. Although this method is common, more accurate results can be achieved by summing up the height of water and the earth subsidence in each scenario and performing probabilistic hazard calculations one time. Finally, we should mention the importance role of high-quality bathymetry and topography data for more accurate PTHA results, especially for simulation of inundation and run-up that we hope to consider in our future studies.
Funding The present study was carried out within the framework of the project "Tsunami and Seismic hazard assessment for the Makran region". Funding for this project was provided by the Plan and Budget Organization of the Islamic Republic of Iran (PBO).