Comprehensive Probabilistic Tsunami Hazard Assessment in the Makran Subduction Zone
CComprehensive Probabilistic Tsunami Hazard Assessment inthe Makran Subduction Zone
Parastoo Salah a , , Jun Sasaki, Mohsen Soltanpour Graduate Program in Sustainability Science-Global Leadership Initiative, Graduate School ofFrontier Sciences, The University of Tokyo, Kashiwa, Chiba, 277-8563, Japan, Department of Socio-Cultural Environmental Studies, Graduate School of Frontier Sciences,The University of Tokyo, Kashiwa, Chiba, 277-8563, Japan Department of Civil Engineering, K. N. Toosi University of Technology, No. 1346, Vali-Asr St.,Tehran, Iran
Abstract
After the 2004 and 2011 tsunamis came unprecedented to the scientific communitythe role of probabilistic tsunami hazard assessment (PTHA) in tsunami-prone areas came to thefore. The Makran subduction zone (MSZ) is a hazardous tsunami-prone region; however, due toits low population density, it is not as prominent in literature. In this study, we assess the threatof tsunami hazard posed to the coast of Iran and Pakistan by the MSZ and present a compre-hensive PTHA for the entire coast regardless of population density. We accounted for sourcesof epistemic uncertainties by employing event tree and ensemble modeling. Aleatory variabilitywas also considered through probability density function. Further, we considered the contribu-tion of small to large magnitudes and used our event trees to create a multitude of scenarios asinitial conditions.
Funwave-TVD was employed to propagate these scenarios. Our results demon-strate that the spread of hazard curves for different locations on the coast is remarkably large,and the probability that a maximum wave will exceed 3 m somewhere along the coast reaches { , , , , } % for return periods { , , , , } , respectively. Moreover, we foundthat the exceedance probability could be higher at the west part of Makran for a long returnperiod, if we consider it as active as the east part of the MSZ. Finally, we demonstrated that thecontribution of aleatory variability is significant, and overlooking it leads to a significant hazardunderestimation, particularly for a long return period. a [email protected] a r X i v : . [ phy s i c s . g e o - ph ] O c t ontents M max : node 3 in EV1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82.2.4 Earthquake catalogues . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92.3 Tsunami scenarios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102.3.1 Source discretization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102.3.2 Rupture area: node 1 and 2 in EV2 . . . . . . . . . . . . . . . . . . . . . . 102.3.3 Slip distribution: node 3 in EV2 . . . . . . . . . . . . . . . . . . . . . . . . 102.3.4 Possible location: node 4 in EV2 . . . . . . . . . . . . . . . . . . . . . . . . 122.4 Tsunami model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122.5 Deriving the probability of exceedance . . . . . . . . . . . . . . . . . . . . . . . . . 122.5.1 Ensemble model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 Tsunami events are infrequent in several water bodies around the world, yet their danger cannotbe ignored due to the high levels of destruction that follow, including major losses of life andproperty damage. In particular, the importance of a comprehensive tsunami hazard assessment(THA) is highlighted when a disastrous tsunami occurs. Recent devastating tsunamis such as theSumatra tsunami of 2004, with more than 200,000 fatalities [1], and the 2011 Tohoku tsunami inJapan, which caused more than 15,000 fatalities and was responsible for the Fukushima NuclearPower Plant accident [2], are representative examples. Following these disasters, there has beena remarkable development in tsunami risk management to reduce the effect of future tsunamis.For recent reviews of these developments, including full lists of references, see [3, 4].Tsunami hazard assessment includes sensitivity analyses (see e.g. [5, 6]) as well as determinis-tic (see e.g. [7, 8, 9]) and probabilistic approaches. The latter approach– called the probabilistic1sunami hazard assessment (PTHA) –has received substantially increased attention after the 2004and 2011 tsunamis [10, 11, 12, 13]. Unlike deterministic approaches that consider specific sce-narios (commonly including the worst case scenario) to calculate tsunami hazard metrics (suchas run up height and arrival time), PTHA calculates the likelihood of tsunami impact employ-ing multiple possible scenarios consisting of the contributions from small to large events alongwith all quantifiable uncertainties [14]. Hence, PTHA can overcome the limitation of incompleteor insufficient historical records, and extend the return periods from hundreds to thousands ofyears. Furthermore, this approach considers the uncertainties stemming from the lack of re-searcher knowledge and the random nature of hazards. The former is represented by epistemicuncertainty in literature while the latter is known as aleatory variability. These concepts areexplained in detail in section 2.1.PTHA was developed by adopting the probabilistic seismic hazard analysis (PSHA) [15, 16,17], and much progress has been built upon it see [18] and the references therein. Notwithstand-ing that PTHA is a relatively new method, it has been widely used in tsunami-prone areas owingto its diverse range of applications (e.g., [19, 20, 21]), each of them covers different uncertain-ties, methods, and level of accuracy. In this study, we assessed the tsunami hazard using theprobabilistic approach in the Makran subduction zone (MSZ).The MSZ is a tsunami risk zones as attested by compiled tsunami catalogues and recentpaleotsunami studies [22] that exhibits risks for the neighboring countries of Iran, Oman, andPakistan. This region is not as prominent in scientific literature as other tsunami-prone subduc-tion zones owing to its low population density, and it remains as one of the least studied regions.The authors of [23] performed the first generation of PTHA in the MSZ. Their results are notreliable for return period far from the typical recurrence time of magnitude M w = 8 . Our methodology aims to calculate the probability of exceeding a /set of tsunami heights atthe Makran coast, considering both epistemic and aleatory uncertainties. In this work, we onlyfocused on tsunamis induced by earthquakes; landslide-induced tsunamis were beyond the scopeof our research and should be addressed in future work. Fig 1 demonstrates a summary of ourframework. First, we determined the seismicity area and generated synthetic scenarios similarto that described by [20]. Then, for each scenario, we ran a fully nonlinear tsunami model
Funwave-TVD [27, 28] to obtain the maximum wave heights along the coastline. Additionally, weincorporated the epistemic uncertainties by developing two event trees and ensemble modeling.Finally, we calculated the tsunami height exceedance rate considering the aleatory variability.
A reliable PTHA must consider the epistemic uncertainty and aleatory variability simultaneously.The latter expresses the innate variability of the physical process, while the former is related tothe lack of understanding and limited knowledge of the process. As shown in Fig 1, each factorwas considered as described in detail below.
Epistemic uncertainties can be incorporated by developing event trees [29]. We developed twoevent trees:(i) Focusing on the fault source recurrence model for the assessment of mean annual rates ofearthquakes at different magnitude levels with 36 branches. It consists of two zonations [seg-mented and none]; three approaches for the seismicity model [Gutenberg–Richter–Bayes (GRB)[30], truncated Gutenberg–Richter, and characteristic [31]]; three maximum magnitudes ( M max )[based on the Kijko–Sellevoll–Bayes method [32], thermomechanical modeling [33] and ergodicassumption [34]]; and three for incorporating the uncertainty of the earthquake occurrence model.See Fig 2(ii) Focusing on the bulk rupture parameters and rupture complexity. It consists of rupturelength and width, earthquake source location within the fault, and slip distribution. See Fig 2 Many authors believe that no theoretical significance exists for this separation because, as long as our knowledgeincreases, all uncertainties become epistemic [26]. igure 1: Methodology framework. First, the fault geometry was defined using SLAB 2.0, and the source wasdiscretized into smaller segments. Next, two event trees were developed to define the earthquake recurrence rateand create tsunami scenarios; then, the Okada model and
Funwave-TVD were used to calculate tsunami heights forour scenarios. Finally, considering the aleatory variability, we derived the probability of exceedance.
The proper treatment of aleatory variability in tsunami wave heights is a prominent subject, andignoring this typically leads to significant hazard underestimation [35]. In our analysis, we haveidentified three main contributions, i.e., { σ m , σ s , σ t } , to the aleatory variability as below. Numerical model and bathymetry ( σ m ) – Due to the lack of field data and backgroundinformation on the MSZ, the 2011 Tohoku earthquake of Japan was modeled, and the resultswere compared with the available measured data to quantify the mismatch between the observedand computed tsunami heights. This uncertainty is described as the standard deviation of a4 a)(b)
Figure 2:
Developed event trees for (a) source recurrence model; (b) rupture complexity and tsunami scenariocreation. log-normal distribution with a zero mean [36, 37]: σ m = log κ = vuut n n X i =1 (log K i ) − (log K ) , log K = 1 n n X i =1 log (cid:18) H obs H model (cid:19) . (1)Here, K i = H obs /H model with H obs and H model are the measured and simulated tsunami heights,respectively. For H obs , the measured tsunami height at GPS, DART buoys, and tide and wavegauges were used. Moreover, we simulated the 2011 Tohoku tsunami using the same bathymetryand numerical model as the ones we used for the MSZ to obtain H model (see section 2.4). Fig 3shows the comparison between the modeled and measured tsunami heights with σ m = 0 . Scaling relations ( σ s ) – Given the earthquake magnitude, rupture length and width werederived by evaluating the scaling relations. To do so, we used Strasser relations [38] as explainedin section 2.3.2. To account for stochasticity in the earthquake dimensions imposed by the scaling The data from these source were used to avoid uncertainties when using survey measuring methods. σ s = 0 . σ s = 0 .
180 for length and width, respectively. The variability in scaling relations was derivedfrom a regression analysis of the relations [38, Table 1]. log ( K i ) ' H Q V L W \ 0 R G H O H G K H L J K W V P 2 E V H U Y H G K H L J K W V P Figure 3:
Comparison between modeled and measured tsunami height for the 2011 Japan tsunami at 15 stationsrecorded by GPS, DART buoys, tide and wave gauges; regression line for modeled versus measured height (bottomright); histogram of errors in log tsunami height and corresponding normal distribution (bottom left).
Tide ( σ t ) – Because the tide level at tsunami arrival time is unknown, tidal variation variabilitymust be included in the PTHA. In the Makran region, the tidal variation is notable, and thepeak-to-peak tidal amplitude is as high as 2-3 m. For this task, we calculated the probability ofexceedance of mean sea level (MSL) from the tidal record at each point of interest (PoI).To calculate tidal record probability, we used a relatively long time-series of record measuredby tidal gauges for each PoI. For PoIs in which a tidal record is not available, we used a linearinterpolation of the closest tidal gauges. This choice seems reasonable because the differences intidal levels along the Makran coast are not significant [39]. Fig 4 illustrates an example of ourmethodology for one PoI, Beris. 6 - - - - - - DateTime1.51.00.50.00.51.01.5 t i d a l s t a g e ( m ) W L G D O V W D J H P S G I Figure 4:
Tidal time series record of one year starting from 2016 for Beris (top); corresponding normal distribution(bottom). For this PoI, σ t = 0 . . The MSZ is located on the southeastern coasts of Iran and southern coasts of Pakistan. Thiszone extends east from the Strait of Hormoz to the Ornch–Nal Fault in Pakistan. It experiencedthe deadliest tsunami that has occurred in the Indian Ocean prior to 2004, and recent smallerearthquakes suggest seismicity on the megathrust. However, poor historical records have led tosignificant uncertainty and complicated hazard potential estimation. Therefore, as mentionedin section 2.1.1, to incorporate the uncertainties associated with the fault source, we developedan event tree (EV1) to assess the mean annual rates ( ν j ) of earthquakes at different magnitudelevels, as described below. The eastern and western parts of the MSZ exhibit extremely different seismicity patterns [40].This, along with its unrecognized bathymetric trench, makes the MSZ a unique subject of analysis.[41] argued that the eastern MSZ is underlain by an oceanic lithosphere, while the western partis possibly underlain by a continental or very low velocity oceanic lithosphere. This, alongwith the more historical seismicity activity at the eastern part, form the hypothesis of east-westsegmentation of the MSZ. However, it remains a controversial issue whether the MSZ should beconsidered segmented in hazard studies because the existence of late Holocene marine terracesalong the eastern and western halves suggests that both can generate megathrust earthquakes [42].We represent both the segmented and non-segmented zone in node 1 of EV1. However, owingto the above mentioned related controversy, we neglected the hypothesis of the segmented MSZ7s it leads to a strong hazard underestimation. Accordingly, we weighted the segmented andnon-segmented zones as 0 and 1, respectively.
The severity of a large earthquake is determined by the tail of a frequency distribution. Thus,earthquake catalogues are limited at large magnitudes for a particular fault zone. This makesthe accurate estimation of the probabilistic tsunami hazard through the application of the re-currence interval of seismic history impossible. In particular, for MSZ with poor and incompletecatalogues, a simple linear regression of the historical cumulative distribution is known to bebiased [43]. Accordingly, several models exist that can be used to define the distribution of earth-quake magnitudes for incomplete catalogues. In this study, we used three seismicity models:(i) Gutenberg–Richter–Bayes(GRB) [30]. Seismicity was determined using the
HA3 applicationbuilt in
MATLAB . The applied procedure of the seismic hazard considers the incompletenessof the seismic catalogues, uncertainty in magnitude estimation, and variation in seismic-ity. The code accepted mixed data catalogues, namely, paleo, historical, and instrumentalwith different completeness magnitudes, time periods, and magnitude uncertainties. Thismethod employed a mixed (Bayesian) Poisson-gamma distribution as a model of earthquakeoccurrence over time.(ii) Characteristic [31]. The characteristic distribution has the cumulative complementary func-tion (Φ M ) truncated on both ends and is characterized by the following equationΦ M = ( e − β ( M − M min ) for M min ≤ M ≤ M max M > M max . (2)(iii) Truncated Gutenberg-Richter (TGR). The cumulative complementary function (Φ M ), whichis truncated at both ends, is expressed asΦ M = ( e − β ( M − M min) − e − β ( M max − M min) − e − β ( M max − M min) for M min ≤ M ≤ M max M > M max , (3)where M min is the level of magnitude completeness, M max is the maximum possible earth-quake magnitude and β = b log 10, and b is the parameter of the Gutenberg-Richter relation. M max : node 3 in EV1 PTHAs are more sensitive to M max than PSHAs because tsunami heights do not saturate withincreasing magnitude as seismic ground motions do [37]. M max based on instrumental cataloguesmay underestimate the maximum magnitude event due to their short records. Here, to includethis uncertainty, we used three methods for maximum magnitude ( M max ) assessment: Note that treating Tohoku as a segmented zone led to strong underestimation of the devastating 2011 tsunami[11]. able 1: Extracted values for magnitude of completeness ( M c ), error, and b -value from zmap for different workingcatalogues. prehistorical historical complete 1 complete 2 complete 3period 326 BC − − − − − M c — 5 . . . . . . .
45 0 .
35 0 . b -value 0 . ± . HA3 application, we found M max = 8 . M max = 9 .
22 for the full length ofsubduction zone in [33].(iii) Ergodic assumption: [34] suggested M max = 9 .
58 for subductions based on their statisticalanalysis for a number of faults worldwide.
Earthquake data employed in this study were derived from various sources: (i) InternationalSeismological Centre (ISC) (ii) Incorporated Research Institutions for Seismology (IRIS) (iii) TheUnited States Geological Survey Online bulletin (USGS), which includes information from theNational Ocean and Atmospheric Administration (NOAA) and Preliminary Determination ofEpicentres (PDE) provided by the National Earthquake Information Center (NEIC) (iv) GlobalHistorical Earthquake Archive (GEM) (v) Iranian Seismological Center (IRSC). Extra effort hasbeen made to extract additional data from literature regarding earthquakes with magnitudesbeyond 6.5. This includes information from the Pakistan Meteorological Department (PMD) [44]and [45].We compiled the catalogues for a region that lies in the plate interface, excluding nonsub-duction seismicity (see Fig. 5). The catalogues cover the period from 825 BCE to mid-2020 CE.These catalogues are different in terms of magnitude scale. When available, the moment mag-nitude, M w , was used; otherwise, the published magnitudes (e.g., teleseismic magnitudes andmodified Mercalli intensity) were converted to M w using the empirical laws proposed by [46, 47].We used the ZMAP7 analysis tool [48] to prepare the catalogues for our recurrence models.First, following the assumption that seismicity obeys a Poisson process, it is necessary to declusterthe catalogues by removing all dependent events, namely, precursors and aftershocks. Hence, weemployed the cluster approach proposed by Reasenberg [49] to eliminate dependent shocks. Then,duplicate events from different catalogues were removed. Sunsequently, the plot of the cumulativenumber of events allowed us to split the working catalogues into prehistorical, historical, and threesub-instrumental categories. Each has a different magnitude of completeness ( M c ) and magnitudeuncertainty. Moreover, we obtained a prior value for b in each catalogue to use in our recurrencemodels (see Table 1). 9 .3 Tsunami scenarios To create possible tsunami scenarios and incorporate rupture and location uncertainties, eventtree 2 (EV2) was developed (see Fig 2). The branches of EV2 are introduced in section 2.1.1;here, we describe them in detail.
Similar to [20], fault geometry was defined using a three-dimensional source zone fault-plane,SLAB 2.0 – a comprehensive subduction zone geometry model [50]. The MSZ has an extremelyshallow subduction angle (dip) and thick sediment pile ( ≈ ×
50 km segments. Finally, dip, rake, strike, and depth for each segmentwere identified for use in the Okada model [53] to generate the initial tsunami conditions. For each magnitude ranging from M w = 7 . to M w = 9 . M w ∈ { M w, min , M w, min + 0 . , . . . , M w, max − . , M w, max } , we calculated the rupturelength and width using the scaling relation of [38, Table 1] derived from the regression analysis ofhistorical subduction events. For M w ≤ .
7, we included uncertainties associated with the use ofthe scaling relation for earthquake dimensions as described in section 2.1. However, we observedthat the variability enlarges with growing magnitude. Hence, for M w > .
7, rather than usingonly one value for rupture length and width, a random sample was selected from a log-normaldistribution. Our initial intention was to consider the dependence of the variance of rupture lengthand width, and sample from a two-dimensional multivariate normal distribution [54]. However,we noticed that the range of variation was quite small. Considering our segmentation size (i.e.,50 ×
50 km ), the independent random selection of length ( L ) and width ( W ) were generatedform normal distributions according tolog L ∼ N ( − .
477 + 0 . M w , . , log W ∼ N ( − .
882 + 0 . M w , . . (4)Here, N ( µ, σ ) is a normally distributed random variable with mean µ and standard deviation σ ;notation ∼ denotes the equivalence of distributions.Then, for each length and width we calculated the number of segments downdip ( n s ) andalong-strike ( n l ) using the method described in [55, Eq. (4)]. Slip distribution significantly affects tsunami heights nearshore. Recently, different studies haveshown that maximum nearshore wave height varies by a factor of 2 or more due to heterogeneity M w = 7 . igure 5: Bathymetry model and computational domain; dots represent PoIs located at 0 m isobath along thecoast; blue mesh indicates the seismogenic zone and source discretization into × km . in earthquake slip [56, 57, 5, 58, 59]. However, owing to its convoluted nature and computationcomplexity, tsunami hazard assessments are usually based on idealized uniform slip earthquakes.In this work, we used a uniform slip for M w ≤ .
7, where the effect of spatial slip distributionis not significant, and heterogeneous slip distribution for M w ≥ .
8, where we found that theheterogeneity of slip notably varies tsunami heights at our PoIs. This trade-off was specifiedto account for the effect on tsunami heights and optimize the number of scenarios through oursensitivity analysis. The number of scenarios and differences among modeled tsunami heights atPoIs were compared for a fixed scenario, but with varying M w , starting from M w = 7 . M w > . M w employing the scaling rela-tion as follows: M w = log M o − . . , S = M o µ × A , (5)where M o is the seismic moment, µ the shear modulus, and A the area of each scenario. We set µ = 3 × dyn cm − as it is appropriate for crustal rocks and shallow depth faults [60]. Wethen used the evaluated S from Eq. (5) as a uniform slip for M w ≤ .
7; whereas for M w > . L , W )-scenario was created randomly using the PTHA18 code built in R . The PTHA18 code uses the SNCF model of [61] for generating random slip distribution for agiven segment dimension and number. This model is a variant on the widely used K model.Further implementation details can be found in [62] and [61].11 .3.4 Possible location: node 4 in EV2 To cover all the seismogenic zone for each magnitude and sampled length and width, we floatedthe calculated n s × n l through all possible locations of the Makran seismicity area, shown inthe blue mesh in Fig. 5. We assumed that the occurrence of a specific magnitude was equallyprobable in all possible locations; therefore, an equal weight was assigned to the branches of EV2. In total, 4220 scenarios were created using the approach discussed in the previous section con-sidering the branches of EV2 for different magnitudes, which were randomly sampled from therupture area, slip distribution, and all possible locations. For each scenario, numerical simulationsof tsunami generation and propagation were performed. First, we calculated vertical co-seismicdislocation via a homogeneous elastic half-space model [53]. Then, the Kajiura filter [63] wasused for the ocean surface deformation of the dislocation to calculate the initial conditions. Re-garding the simulation of tsunami propagation a fully nonlinear and dispersive Boussinesq longwave model,
FUNWAVE-TVD [27, 28], was employed. It features accurate dissipation by consideringthe breaking wave and bottom friction processes, and has been systematically validated againstexperimental studies and benchmarks [64]. The code was parallelized using the message pass-ing interface (
MPI ). This salable algorithm (using more than 90% of the number of cores in acomputer cluster [28]) has been paved our way for modeling multitude scenarios.Here,
FUNWAVE-TVD was used in its Cartesian implementation. To prevent non-physical reflec-tion from the boundaries, sponge layers were specified with 10 km thickness within our computa-tional domain. A 600 m resolution was used for the computational domain, which is a trade-offbetween precision and practical feasibility.To guarantee representative bathymetry, we evaluated three bathymetric models. In par-ticular, Etopo-v1 (based on satellite gravity data), GEBCO (global bathymetric model basedon ship-track data), and SRTM+ (space shuttle radar mapping) with measured data providedby the Ports and Maritime Organization of Iran (PMO) were compared. It appears that thelatest released data of the GEBCO model with 15 arcsecond resolution are the best among theaforementioned models and exhibit a smoother transition between deep and shallow water. Thiscontrasts the results of [37] while agreeing with those of [65]. The former can be due to the het-erogeneous characteristic of each site or/and the latest update of the bathymetry data. Hence,GEBCO-2020 has been used for our tsunami simulations. Each scenario has been simulated for8 h, and for each computational time step, a time series of tsunami wave has been recorded at 84hazard points. These PoIs are located at 5 to 0 m isobath at approximately 10 −
12 km intervalsalong the Iran and Pakistan coastline. The PoIs are shown in Fig 5.
For a given exposure time (∆ T ), PTHA was performed by deriving the exceedance of maximumtsunami height ( ψ ) at each PoI form a threshold value ( ψ t ). Considering a total of J possible12agnitudes, we defined the total probability of exceedance P tot ( ψ > ψ t , ∆ T, PoI) = 1 − J Y j =1 (1 − P ( E j , ∆ T ) P ( ψ > ψ t | E j )) , (6)where P ( E j , ∆ T ) is the probability that at least one event ( E j ) occurs in the return period ∆ T .Assuming that the occurrence of earthquakes conforms to a stationary Poisson process with theannual recurrence rate ν j , it can be assessed as P ( E j , ∆ T ) = 1 − exp ( − ν j × ∆ T ) . (7)Considering the uncertainties on rupture dimensions, locations, and slip distribution in EV2,each E j can cause different scenarios ( S ( j ) A ). The probability that tsunami height ( ψ ) exceeds athreshold ( ψ t ) when the event E j occurs is then given by P ( ψ > ψ t | E j ) = A j X A =1 P (cid:0) S ( j ) A | E j (cid:1) P (cid:0) ψ > ψ t |S ( j ) A (cid:1) . (8)Here, P (cid:0) S ( j ) A | E j (cid:1) is the probability of occurrence of the scenario S ( j ) A , and in the absence ofaleatory variability, P (cid:16) ψ > ψ t (cid:12)(cid:12)(cid:12) S ( j ) A (cid:17) = (cid:26) , ψ < ψ t , ψ ≥ ψ t , (9)While in the presence of the aleatory variability that was discussed in section 2.1.2 [66], P (cid:16) ψ > ψ t (cid:12)(cid:12)(cid:12) S ( j ) A (cid:17) = 1 − Φ (cid:18) log( ψ t ) (cid:12)(cid:12) (cid:20) log( ψ S ( j ) A ) (cid:21) , σ (cid:19) . (10)Φ is the cumulative distribution function for a log-normal distribution with the mean equal tothe modeled tsunami height at each PoI and standard deviation σ , given value of a log( ψ t ).From [66], σ can be computed by combining our aleatory variability terms σ = q σ + σ + σ . In this section, we explain how to incorporate the uncertainties from EV1 and obtain P ( E j , ∆ T )using the ensemble model [26]. To calculate the probability that at least one earthquake E j occurs for the selected ∆ T , an event tree was developed as described in section 2.2 and Fig 2.The branches of EV1 were treated in the framework of ensemble modeling, as introduced in [26].Ensemble modeling presumes that epistemic uncertainty is greater than that evaluated by anevent tree, and treats the branches of the event tree as an unbiased sample from a parent distri-bution. This distribution, f ( θ ), describes the variable θ simultaneously considering the aleatoryvariability and epistemic uncertainty.In our case, branches of EV1 are a small sample size, and their few probability outcomescan be replaced by a parametric distribution. A natural choice is the beta distribution that iscommonly used in hazard literature. In this case, we set variable θ ( E j ) = P ( E j , ∆ T ) so that thevariable will be the hazard curve. Different θ ( E j ) are the branches of EV1 that are now a sample13f a Beta ( α, β ) distribution. Parameters α and β are related to the average and variance of θ ( E j ) as E[ θ ( E j ) ] = αα + β , Var[ θ ( E j ) ] = αβ ( α + β ) ( α + β + 1) . (11)In our context, E[ θ ( E j ) ] and Var[ θ ( E j ) ] denote, respectively, the weighted average and variance ofthe exceedance probabilities of the j th magnitude for the selected ∆ T . Inverting equations (11),we found the parameters of the Beta distribution for each magnitude j . Finally, calculating theBeta parameters of the exceedance probability for a set of magnitudes, we plotted the full hazardcurve. In this section, we present the results obtained from the analyses and modeling presented inprevious sections for the coastal area of the MSZ. Our main results are presented by earthquakeand tsunami probability exceedance curves and tsunami probability maps for the selected returntime periods. In this study, we set ∆ T = { , , , , } years; each choice interestsdifferent stakeholders and provides information on a specific aspect of the tsunami hazard inthe MSZ. We also compared the results obtained in the presence and absence of the aleatoryvariability. The earthquake probability of exceedance for the selected ∆ T s are depicted in Fig 6. Theensemble model results from section 2.5.1 is shown through its statistical description, its mean,and the 16th-86th percentiles confidence intervals. For the sake of comparison, the branches ofEV1 ( θ ( E j ) ) are also shown in light gray. In nearly all cases, the statistical description of themean ensemble model is a good representative of EV1 branches Fig 6. Henceforth, we use thevalue of mean ensemble for each magnitude to calculate tsunami probability exceedance curvesand probability maps. Using the equations described in section 2.5, we calculated the probability of exceedance froma set of tsunami height thresholds ψ t = { . , , . , .., . , } m and generated hazard curves atdifferent PoIs incorporating all uncertainties described in the previous sections. The results aregiven in Fig. 7. The hazard curves for each PoI are shown in gray. The results show that thespread of hazard curves for different locations of the Makran coast is remarkably large. As anexample, P tot ( ψ > , ∆ T = 50) ranges from 0 to 16% for different PoIs. By increasing ∆ T , thisrange opens up and it reaches 30% , , M w - y e a r s E x c ee d a n c e P r o b a b ili t y EV1 branchesEnsemble_meanEnsemble_84 th ,16 th percentile M w - y e a r s E x c ee d a n c e P r o b a b ili t y EV1 branchesEnsemble_meanEnsemble_84 th ,16 th percentile M w - y e a r s E x c ee d a n c e P r o b a b ili t y EV1 branchesEnsemble_meanEnsemble_84 th ,16 th percentile M w - y e a r s E x c ee d a n c e P r o b a b ili t y EV1 branchesEnsemble_meanEnsemble_84 th ,16 th percentile M w - y e a r s E x c ee d a n c e P r o b a b ili t y EV1 branchesEnsemble_meanEnsemble_84 th ,16 th percentile Figure 6:
Earthquake probability of exceedance for our sample of ∆ T s: red and blue curves show the statisticaldescription of ensemble model, i.e., mean and th- th percentiles, respectively. For comparison, all outcomes ofthe EV1 branches are also displayed in light gray. period, the probability of exceedance does not vary much among different cities. However, thisdifference becomes significant with increasing ∆ T . For instance, for P tot ( ψ > , ∆ T = 1000), itranges from 32% in the west (Gwadar) to 97% in the east (Jask). Moreover, the probability thattsunami height exceeds 4 m is low (less than 10%) near all major cities except for Jask, with theprobability of 53%. 15 t - y e a r s E x c ee d a n c e P r o b a b ili t y t - y e a r s E x c ee d a n c e P r o b a b ili t y t - y e a r s E x c ee d a n c e P r o b a b ili t y t - y e a r s E x c ee d a n c e P r o b a b ili t y t - y e a r s E x c ee d a n c e P r o b a b ili t y Figure 7:
Tsunami probability of exceedance for a sample of ∆ T s at different PoIs along the Iran and Pakistancoasts. The effect of inclusion of the aleatory variability introduced in section 2.1.2 is shown in Fig. 9.Fig. 9 (a) illustrates the probability of exceedance in the presence and absence of the aleatoryvariability at one random PoI (i.e., Chabahar) for two return periods (100 and 500 years). Theinclusion of the aleatory variability has a significant effect on the probability of exceedance,which increases for a longer return period. As an example, the differences between P tot ( ψ > , ∆ T, Chabahar) with and without the aleatory variability are 8% and 26% for ∆ T = 100 and500 years, respectively. To obtain a better interpretation, we also calculated this difference forall PoIs and ∆ T s, see Fig. 9 (b). In summary, omitting the aleatory variability mostly leads to a16 t E x c ee d a n c e P r o b a b ili t y a t C h a b a h a r
50 years100 years250 years500 years1000 years 0 1 2 3 4 5 6 t E x c ee d a n c e P r o b a b ili t y a t G w a d a r
50 years100 years250 years500 years1000 years0 1 2 3 4 5 6 t E x c ee d a n c e P r o b a b ili t y a t J i w a n i
50 years100 years250 years500 years1000 years 0 1 2 3 4 5 6 t E x c ee d a n c e P r o b a b ili t y a t K o n a r a k
50 years100 years250 years500 years1000 years0 1 2 3 4 5 6 t E x c ee d a n c e P r o b a b ili t y a t J a s k
50 years100 years250 years500 years1000 years 0 1 2 3 4 5 6 t E x c ee d a n c e P r o b a b ili t y a t R a m i n
50 years100 years250 years500 years1000 years
Figure 8:
Tsunami probability of the selected ∆ T s exceedance for the selected six PoIs near main cities. noticeable underestimation with a median of 10% for all PoIs, reaching and it even reaches 40%at somewhere for 1000-year return period. We used a probability map to assess hazard along the entire coast irrespective of populationdensity, which is crucial for prioritizing tsunami mitigation plans and city development in low-population areas. Most literature and mitigation plans focus on specific populated areas. Fig. 10illustrates tsunami probability maps exceeding from two selective thresholds, ψ t = 1 , .0 0.5 1.0 1.5 2.0 2.5 3.0 t - y e a r s E x c ee d a n c e P r o b a b ili t y a t C h a b a h a r without aleatory variabilitywith aleatory variability 0.0 0.5 1.0 1.5 2.0 2.5 3.0 t - y e a r s E x c ee d a n c e P r o b a b ili t y a t C h a b a h a r without aleatory variabilitywith aleatory variability (a) T=50 T=100 T=250 T=500 T=1000Different return periods (years)10010203040 E x c ee d a n c e p r o b a b ili t y d i fff e r e n c e s ( % ) (b) Figure 9: (a) Exceedance curve of Chabahar as a random PoI for 100- and 500-year return periods; blue and redcurves show the probability of exceedance in the presence and absence of the aleatory variability, respectively; (b)box plot showing the differences in exceedance probability ( % ) for different ∆ T s with and without the presence ofthe aleatory variability for all PoIs. tsunami height will exceed 3 m for return periods of 100 and 1000 years is approximately 30%and 95%, respectively. Notably, this is almost 6 to 7 times higher than that in Chabahar. Owingto the small distances between these regions, the inundated area at Chabahar may be affected.Inundation maps are beyond the scope of this study, and we plan to address them in a futurework. The MSZ is one of the two sources of tsunamis in the Indian Ocean, and has the potential of gener-ating large tsunamis that threaten neighboring countries of Iran, Oman, and Pakistan. However,a fortune lack of large tsunamis recently has led to a false sense of safety between community18 La t i t ude Longitude P tot ( >1, T= 1000)
50 mi 100 km La t i t ude Longitude P tot ( >3, T= 1000)
50 mi 100 km La t i t ude Longitude P tot ( >1, T= 500)
50 mi 100 km La t i t ude Longitude P tot ( >3, T= 500)
50 mi 100 km La t i t ude Longitude P tot ( >1, T= 250)
50 mi 100 km La t i t ude Longitude P tot ( >3, T= 250)
50 mi 100 km La t i t ude Longitude P tot ( >1, T= 100)
50 mi 100 km La t i t ude Longitude P tot ( >3, T= 100)
50 mi 100 km
Figure 10:
Maps of tsunami probability exceeding 1 and 3 m for different ∆ T s along the entire coast of Iran andPakistan. T = { , , , , } ranges from 0 to { , , , , } percent, respectively, for differentPoIs.2. The probability of exceedance at PoIs near populated cities decreases and becomes insignif-icant for the exceedance threshold of 4 m (even for a long return period), except for Jaskat the western coast of Iran. Our results provide evidence that if we consider the westernpart of the MSZ equally active and potential to the eastern part –similar to this paper, byweighting both parts the same in our event tree–the exceedance probability could be higherat the western part for a long return period. This can be clearly seen from the probabilitymaps where the exceedance probability of 3 m fluctuates and becomes maximum at thewestern part of the MSZ.3. The inclusion of the aleatory variability has a significant effect on the probability of ex-ceedance, and not including it mostly leads to a remarkable underestimation in the PTHAwith a median of 10% difference for all PoIs. This difference is underscored by increasing thereturn periods and reaches 40% at somewhere for 1000-year return period in the presenceand absence of the aleatory variability.Owing to Makran’s economical, geographical, and strategic importance, Iran approved a planfor developing the southern Makran Coast on December, 2016 titled “Makran Sustainable De-velopment.” This plan, along with the drought occurring recently in the neighboring cities, hasled to an inevitable migration toward the coastlines, with Chabahar exhibiting 10% populationrate growth last year and ranking among the highest population growth rates globally in 2019.Hence, our results are of vital for various stakeholders for developing and implementing tsunamirisk mitigation activities and guiding risk-aware city planning.This study is the first step toward comprehensive and reliable mitigation plans and activities;however, it is important to acknowledge its limitations. Tsunami sources beyond earthquakes were20ot considered in this study. Notably, the only tsunami induced with combination of earthquakeand landslide in word history occurred in Makran in 1945. Moreover, in September 2013, alandslide was recorded immediately following an earthquake in the MSZ [67]. This highlights theneed to consider landslides [68] and their combination with earthquakes in future PTHA studies.Moreover, we disregarded the dynamic interaction between tides and tsunami waves. This worksfor a tsunami wave with one isolated peak; however, it may lead to hazard underestimation whenthe tsunami has several peaks with significant heights. Finally, providing an accurate inundationmap is paramount, which is the aim of our next study. Acknowledgments
References [1] R. Starrs,
When the tsunami came to shore: Culture and disaster in Japan . Global Oriental,2014.[2] Y. Fujii, K. Satake, S. Sakai, M. Shinohara, and T. Kanazawa, “Tsunami source of the 2011off the pacific coast of tohoku earthquake,”
Earth, planets and space , vol. 63, no. 7, p. 55,2011.[3] F. Løvholt, J. Griffin, and M. Salgado-G´alvez, “Tsunami hazard and risk assessment on theglobal scale,”
Encyclopedia of complexity and systems science , pp. 1–34, 2015.[4] P. J. Ward, V. Blauhut, N. Bloemendaal, J. E. Daniell, M. C. de Ruiter, M. J. Duncan,R. Emberson, S. F. Jenkins, D. Kirschbaum, M. Kunz, et al. , “Natural hazard risk assess-ments at the global scale,”
Natural Hazards and Earth System Sciences , vol. 20, no. 4, 2020.[5] K. Goda, P. M. Mai, T. Yasuda, and N. Mori, “Sensitivity of tsunami wave profiles and inun-dation simulations to earthquake slip and fault geometry for the 2011 tohoku earthquake,”
Earth, Planets and Space , vol. 66, no. 1, p. 105, 2014.[6] K. Goda, N. Mori, T. Yasuda, A. Prasetyo, A. Muhammad, and D. Tsujio, “Cascadinggeological hazards and risks of the 2018 sulawesi indonesia earthquake and sensitivity analysisof tsunami inundation simulations,”
Frontiers in Earth Science , vol. 7, p. 261, 2019.[7] P. Lynett, Y. Wei, and D. R. Arcas,
Tsunami Hazard Assessment: Best Modeling Practicesand State-of-the Art Technology . United States Nuclear Regulatory Commission, Office ofNuclear Regulatory . . . , 2016. 218] M. Heidarzadeh, M. D. Pirooz, and N. H. Zaker, “Modeling the near-field effects of theworst-case tsunami in the makran subduction zone,”
Ocean Engineering , vol. 36, no. 5,pp. 368–376, 2009.[9] P. Salah and M. Soltanpour, “Modeling of tsunami propagation and inundation at the northcoast of gulf of oman due to earthquake in makran subduction zone,”
International conferenceon coasts, ports and marine structures (icopmas) , vol. 11, 2014.[10] K. Satake, “Advances in earthquake and tsunami sciences and disaster risk reduction sincethe 2004 indian ocean tsunami,”
Geoscience Letters , vol. 1, no. 1, p. 15, 2014.[11] Y. Y. Kagan and D. D. Jackson, “Tohoku earthquake: A surprise?,”
Bulletin of the Seismo-logical Society of America , vol. 103, no. 2B, pp. 1181–1194, 2013.[12] S. Lorito, J. Selva, R. Basili, F. Romano, M. Tiberti, and A. Piatanesi, “Probabilistic hazardfor seismically induced tsunamis: accuracy and feasibility of inundation maps,”
GeophysicalJournal International , vol. 200, no. 1, pp. 574–588, 2015.[13] F. Løvholt, N. J. Setiadi, J. Birkmann, C. B. Harbitz, C. Bach, N. Fernando, G. Kaiser, andF. Nadim, “Tsunami risk reduction–are we better prepared today than in 2004?,”
Interna-tional journal of disaster risk reduction , vol. 10, pp. 127–142, 2014.[14] E. L. Geist and P. J. Lynett, “Source processes for the probabilistic assessment of tsunamihazards,”
Oceanography , vol. 27, no. 2, pp. 86–93, 2014.[15] C. Cornell, “Engineering seismic risk analysis, bulletin of the seismological society of amer-ica,” 1968.[16] T. Rikitake and I. Aida, “Tsunami hazard probability in japan,”
Bulletin of the SeismologicalSociety of America , vol. 78, no. 3, pp. 1268–1278, 1988.[17] G. L. Downes and M. W. Stirling, “Groundwork for development of a probabilistic tsunamihazard model for new zealand,” in
International Tsunami Symposium 2001 , pp. 293–301,Pacific Marine Environmental Lab. Seattle, Wash., 2001.[18] A. Grezio, A. Babeyko, M. A. Baptista, J. Behrens, A. Costa, G. Davies, E. L. Geist,S. Glimsdal, F. I. Gonz´alez, J. Griffin, et al. , “Probabilistic tsunami hazard analysis: Multiplesources and global applications,”
Reviews of Geophysics , vol. 55, no. 4, pp. 1158–1198, 2017.[19] H. K. Thio, R. I. Wilson, and K. Miller, “Evaluation and application of probabilistic tsunamihazard analysis in california,”
AGUFM , vol. 2014, pp. NH12A–01, 2014.[20] G. Davies and J. Griffin, “Sensitivity of probabilistic tsunami hazard assessment to far-fieldearthquake slip complexity and rigidity depth-dependence: Case study of australia,”
Pureand Applied Geophysics , pp. 1–28, 2019.[21] N. Mori, K. Goda, and D. Cox, “Recent process in probabilistic tsunami hazard analy-sis (ptha) for mega thrust subduction earthquakes,” in
The 2011 Japan Earthquake andTsunami: Reconstruction and Restoration , pp. 469–485, Springer, 2018.2222] D. Kakar, G. Naeem, A. Usman, A. Mengal, A. Naderi Beni, M. Afarin, H. Ghaffari, H. Fritz,F. Pahlevan, E. Okal, et al. , “Remembering the 1945 makran tsunami; interviews withsurvivors beside the arabian sea,”
UNESCO-IOC Brochure , vol. 1, p. 79, 2015.[23] M. Heidarzadeh and A. Kijko, “A probabilistic tsunami hazard assessment for the makransubduction zone at the northwestern indian ocean,”
Natural hazards , vol. 56, no. 3, pp. 577–593, 2011.[24] A. Hoechner, A. Y. Babeyko, and N. Zamora, “Probabilistic tsunami hazard assessment forthe makran region with focus on maximum magnitude assumption,”
Natural Hazards andEarth System Sciences (NHESS) , vol. 16, pp. 1339–1350, 2016.[25] I. El-Hussain, R. Omira, A. Deif, Z. Al-Habsi, G. Al-Rawas, A. Mohamad, K. Al-Jabri, andM. A. Baptista, “Probabilistic tsunami hazard assessment along oman coast from submarineearthquakes in the makran subduction zone,”
Arabian Journal of Geosciences , vol. 9, no. 15,p. 668, 2016.[26] W. Marzocchi, M. Taroni, and J. Selva, “Accounting for epistemic uncertainty in psha: Logictree and ensemble modeling,”
Bulletin of the Seismological Society of America , vol. 105, no. 4,pp. 2151–2159, 2015.[27] J. T. Kirby, G. Wei, Q. Chen, A. B. Kennedy, and R. A. Dalrymple, “Funwave 1.0: fullynonlinear boussinesq wave model-documentation and user’s manual,” research report NO.CACR-98-06 , 1998.[28] F. Shi, J. T. Kirby, J. C. Harris, J. D. Geiman, and S. T. Grilli, “A high-order adaptivetime-stepping tvd solver for boussinesq modeling of breaking waves and coastal inundation,”
Ocean Modelling , vol. 43, pp. 36–51, 2012.[29] T. Annaka, K. Satake, T. Sakakiyama, K. Yanagisawa, and N. Shuto, “Logic-tree approachfor probabilistic tsunami hazard analysis and its applications to the japanese coasts,” in
Tsunami and its hazards in the Indian and Pacific Oceans , pp. 577–592, Springer, 2007.[30] A. Kijko, A. Smit, and M. A. Sellevoll, “Estimation of Earthquake Hazard Parameters fromIncomplete Data Files. Part III. Incorporation of Uncertainty of Earthquake-OccurrenceModel,”
Bulletin of the Seismological Society of America , vol. 106, pp. 1210–1222, 05 2016.[31] Y. Y. Kagan, “Seismic moment distribution revisited: I. statistical results,”
GeophysicalJournal International , vol. 148, no. 3, pp. 520–541, 2002.[32] A. Kijko, “Estimation of the maximum earthquake magnitude, m max,”
Pure and AppliedGeophysics , vol. 161, no. 8, pp. 1655–1681, 2004.[33] G. L. Smith, L. C. McNeill, K. Wang, J. He, and T. J. Henstock, “Thermal structure andmegathrust seismogenic potential of the makran subduction zone,”
Geophysical ResearchLetters , vol. 40, no. 8, pp. 1528–1533, 2013.2334] P. Bird and Y. Y. Kagan, “Plate-tectonic analysis of shallow seismicity: Apparent boundarywidth, beta, corner magnitude, coupled lithosphere thickness, and coupling in seven tectonicsettings,”
Bulletin of the Seismological Society of America , vol. 94, no. 6, pp. 2380–2399,2004.[35] J. J. Bommer and N. A. Abrahamson, “Why do modern probabilistic seismic-hazard analysesoften lead to increased hazard estimates?,”
Bulletin of the Seismological Society of America ,vol. 96, no. 6, pp. 1967–1977, 2006.[36] I. Aida, “Reliability of a tsunami source model derived from fault parameters,”
Journal ofPhysics of the Earth , vol. 26, no. 1, pp. 57–73, 1978.[37] H. K. Thio, P. Somerville, and G. Ichinose, “Probabilistic analysis of strong ground motionand tsunami hazards in southeast asia,”
Journal of Earthquake and Tsunami , vol. 1, no. 02,pp. 119–137, 2007.[38] F. O. Strasser, M. Arango, and J. J. Bommer, “Scaling of the source dimensions of interfaceand intraslab subduction-zone earthquakes with moment magnitude,”
Seismological ResearchLetters , vol. 81, no. 6, pp. 941–950, 2010.[39] P. Akbari, M. Sadrinasab, V. Chegini, and S. Siadat Mousavi, “Study of tidal componentsamplitude distribution in the persian gulf, gulf of oman and arabian sea using numericalsimulation,”
Journal of Marine Science and Technology , vol. 16, no. 3, pp. 27–41, 2017.[40] G. Aldama-Bustos, J. Bommer, C. Fenton, and P. Stafford, “Probabilistic seismic hazardanalysis for rock sites in the cities of abu dhabi, dubai and ra’s al khaymah, united arabemirates,”
Georisk , vol. 3, no. 1, pp. 1–29, 2009.[41] A. I. Al-Lazki, K. S. Al-Damegh, S. Y. El-Hadidy, A. Ghods, and M. Tatar, “Pn-velocitystructure beneath arabia–eurasia zagros collision and makran subduction zones,”
GeologicalSociety, London, Special Publications , vol. 392, no. 1, pp. 45–60, 2014.[42] R. Normand, G. Simpson, F. Herman, R. H. Biswas, and A. Bahroudi, “Holocene sedimen-tary record and coastal evolution in the makran subduction zone (iran),”
Quaternary , vol. 2,no. 2, p. 21, 2019.[43] W. Power, L. Wallace, X. Wang, and M. Reyners, “Tsunami hazard posed to new zealandby the kermadec and southern new hebrides subduction margins: an assessment based onplate boundary kinematics, interseismic coupling, and historical seismicity,”
Pure and appliedgeophysics , vol. 169, no. 1-2, pp. 1–36, 2012.[44] P. P. M. Department], “Seismic hazard analysis and zonation for pakistan, azad jammu andkashmir,” 2007.[45] N. Ambraseys and C. Melville, “A history of persian earthquakes cambridge univ,”
Press,New York , 1982. 2446] T. I. Allen, D. J. Wald, and C. B. Worden, “Intensity attenuation for active crustal regions,”
Journal of seismology , vol. 16, no. 3, pp. 409–433, 2012.[47] B. Lolli, P. Gasperini, and G. Vannucci, “Empirical conversion between teleseismic magni-tudes (mb and m s) and moment magnitude (m w) at the global, euro-mediterranean anditalian scale,”
Geophysical Journal International , vol. 199, no. 2, pp. 805–828, 2014.[48] C. Reyes and S. Wiemer, “Zmap7, a refreshed software package to analyze seismicity.,” in
Geophysical Research Abstracts , vol. 21, 2019.[49] P. Reasenberg, “Second-order moment of central california seismicity, 1969–1982,”
Journalof Geophysical Research: Solid Earth , vol. 90, no. B7, pp. 5479–5495, 1985.[50] G. P. Hayes, G. L. Moore, D. E. Portner, M. Hearne, H. Flamme, M. Furtney, and G. M.Smoczyk, “Slab2, a comprehensive subduction zone geometry model,”
Science , vol. 362,no. 6410, pp. 58–61, 2018.[51] K. Berryman, L. Wallace, G. Hayes, P. Bird, K. Wang, R. Basili, T. Lay, M. Pagani, R. Stein,T. Sagiya, et al. , “The gem faulted earth subduction interface characterisation project, ver-sion 2.0, april 2015,”
Global Earthquake Model , 2015.[52] A. Safari, A. Abolghasem, N. Abedini, and Z. Mousavi, “Assessment of optimum value fordip angle and locking rate parameters in makran subduction zone.,”
International Archivesof the Photogrammetry, Remote Sensing & Spatial Information Sciences , vol. 42, 2017.[53] Y. Okada, “Surface deformation due to shear and tensile faults in a half-space,”
Bulletin ofthe seismological society of America , vol. 75, no. 4, pp. 1135–1154, 1985.[54] L. Blaser, F. Kr¨uger, M. Ohrnberger, and F. Scherbaum, “Scaling relations of earthquakesource parameter estimates with special focus on subduction environment,”
Bulletin of theSeismological Society of America , vol. 100, no. 6, pp. 2914–2926, 2010.[55] G. Davies, J. Griffin, F. Løvholt, S. Glimsdal, C. Harbitz, H. K. Thio, S. Lorito, R. Basili,J. Selva, E. Geist, et al. , “A global probabilistic tsunami hazard assessment from earthquakesources,”
Geological Society, London, Special Publications , vol. 456, no. 1, pp. 219–244, 2018.[56] F. Løvholt, G. Pedersen, S. Bazin, D. K¨uhn, R. E. Bredesen, and C. Harbitz, “Stochasticanalysis of tsunami runup due to heterogeneous coseismic slip and dispersion,”
Journal ofGeophysical Research: Oceans , vol. 117, no. C3, 2012.[57] C. Mueller, W. Power, S. Fraser, and X. Wang, “Effects of rupture complexity on localtsunami inundation: Implications for probabilistic tsunami hazard assessment by example,”
Journal of Geophysical Research: Solid Earth , vol. 120, no. 1, pp. 488–502, 2015.[58] H. Sugino, Y. Iwabuchi, N. Hashimoto, K. Matsusue, K. Ebisawa, H. Kameda, and F. Ima-mura, “The characterizing model for tsunami source regarding the inter-plate earthquaketsunami,”
Journal of JAEE , vol. 15, no. 3, pp. 3 114–3 133, 2015.2559] R. Butler, D. Walsh, and K. Richards, “Extreme tsunami inundation in hawai ‘i fromaleutian–alaska subduction zone earthquakes,”
Natural Hazards , vol. 85, no. 3, pp. 1591–1619, 2017.[60] A. Deif and I. El-Hussain, “Seismic moment rate and earthquake mean recurrence intervalin the major tectonic boundaries around oman,”
Journal of Geophysics and Engineering ,vol. 9, no. 6, pp. 773–783, 2012.[61] G. Davies, N. Horspool, and V. Miller, “Tsunami inundation from heterogeneous earthquakeslip distributions: Evaluation of synthetic source models,”
Journal of Geophysical Research:Solid Earth , vol. 120, no. 9, pp. 6431–6451, 2015.[62] G. Davies, “Tsunami variability from uncalibrated stochastic earthquake models: testsagainst deep ocean observations 2006–2016,”
Geophysical Journal International , vol. 218,no. 3, pp. 1939–1960, 2019.[63] K. Kajiura, “The leading wave of a tsunami,”
Bulletin of the Earthquake Research Institute,University of Tokyo , vol. 41, no. 3, pp. 535–571, 1963.[64] B. Tehranirad, F. Shi, J. T. Kirby, J. C. Harris, and S. Grilli, “Tsunami benchmark resultsfor fully nonlinear boussinesq wave model funwave-tvd, version 1.0,”
Center for AppliedCoastal Research, University of Delaware, Tech. Rep , 2011.[65] K. Marks and W. Smith, “An evaluation of publicly available global bathymetry grids,”
Marine Geophysical Researches , vol. 27, no. 1, pp. 19–34, 2006.[66] H. K. Thio, “Urs probabilistic tsunami hazard system: A user manual,”
Technical ReportURS Corporation, San Francisco, CA. , 2010.[67] M. Heidarzadeh and K. Satake, “Possible sources of the tsunami observed in the northwest-ern indian ocean following the 2013 september 24 m w 7.7 pakistan inland earthquake,”
Geophysical Journal International , vol. 199, no. 2, pp. 752–766, 2014.[68] E. Rastgoftar and M. Soltanpour, “Study and numerical modeling of 1945 makran tsunamidue to a probable submarine landslide,”