Probabilistic tsunami hazard assessment for the Makran region with focus on maximum magnitude assumption

Despite having been rather seismically quiescent for the last decades, the Makran subduction zone is capable of hosting destructive earthquakes and tsunami. In particular, the well-known thrust event in 1945 (Balochistan earthquake) led to about 4000 casualties. Nowadays, the coastal regions are more densely populated and vulnerable to similar events. Furthermore, some recent publications discuss rare but significantly larger events at the Makran subduction zone as possible scenarios. We analyze the instrumental and historical seismicity at the subduction plate interface and generate various synthetic earthquake catalogs spanning 300 000 years with varying magnitude-frequency relations. For every event in the catalogs we compute estimated tsunami heights and present the resulting tsunami hazard along the coasts of Pakistan, Iran and Oman in the form of probabilistic tsunami hazard curves. We show how the hazard results depend on variation of the Gutenberg–Richter parameters and especially maximum magnitude assumption.


Introduction
A subduction zone along the Makran coast was proposed for the first time 40 years ago (Stoneley, 1974) and found to be in agreement with seismicity (Quittmeyer and Jacob, 1979) as well field observations (Page et al., 1979).There is no clear trench developed where the Arabian plate subducts beneath the Eurasian plate, and rather low seismicity characterizes the interplate seismogenic zone.Its largest recorded earthquake was an event with M W = 8.1 in 1945 (Byrne et al., 1992).For these reasons, and due to its low population density, it is not as prominent in scientific literature as other tsunami-prone subduction zones.Nevertheless, the Makran or Balochistan event mentioned above caused the tsunami with the highest number of fatalities of about 4000 in the Indian Ocean region prior to the 2004 Sumatra-Andaman earthquake (Heidarzadeh et al., 2008b).
Some publications assessed the tsunami hazard for the coast of Iran, Pakistan and Oman originating from the subduction process: Heidarzadeh et al. (2008a, c) performed deterministic analyses based on five magnitude 8.1 events and six magnitude 8.3 events, respectively, and Heidarzadeh and Kijko (2010) carried out a probabilistic study based on three magnitude 8.1 sources associated with probabilities to propose what they labeled a "first generation" probabilistic tsunami hazard assessment (PTHA) for the region.
The two largest earthquakes of the past 50 years, the 2004 Sumatra-Andaman earthquake, and the 2011 Tohoku earthquake both defied the expectations of most of the scientific community: since such huge events had not been observed before at these locations, they had not been accounted for.Since then, increased attention has been given to the question of what maximum magnitude can be expected for the different subduction zones.For instance, Schellart and Rawlinson (2013) tried to find correlations between 24 parameters related to subduction (geometry, geology, kinematics, dynamics) and maximum magnitude on a global scale, and found that some parameters appear to lie in a limited band for subduction segments capable of very large events.They concluded that Makran falls into the set of regions that could generate events larger than M W = 8.5.McCaffrey (2008) states that present evidence cannot rule out that any subduction zone may produce a magnitude 9 event, based on statistical analysis of Monte Carlo simulations.Smith et al. (2013) performed a thermomechanical analysis to assess the seismogenic potential of the Makran subduction zone.They found that the thick sediment cover leading to high plate boundary temperatures at the deformation front makes the megathrust potentially seismogenic to a shallow depth, and the low dip angle leads to a wide seismogenic zone, resulting in an upper estimate for a seaward earthquake of M W 8.7-9.2.
In this study, we perform a PTHA using synthetic earthquake catalogs based on an analysis of the seismicity in order to capture the hazard from the full spectrum of possible earthquakes for the coastlines of Iran, Pakistan and Oman.Since the recent studies mentioned above suggest that maximum magnitude for Makran might be significantly higher than previously thought, special focus is on the consequences of the assumptions made for maximum magnitude.

Methods
Our approach is similar to the one described by Sørensen et al. (2012): first we determine the seismicity in the area of interest; then we generate synthetic earthquake catalogs spanning a much larger time than instrumental and even historic records in order to get stable statistics.For every single event we run a tsunami model to get estimated maximum wave heights along the coastline.These are used to compute tsunami hazard and finally to generate various hazard plots.The difference to the approach mentioned above is, that instead of defining randomly distributed single-fault events of a certain parameter range within the source areas, we project the synthetic earthquake catalog hypocenters to the subduction plate interface and employ a realistic slip distribution.The procedure is described in more detail in the following sections.

Seismicity
From the various published catalogs for the Middle East we select the most recent one by Zare et al. (2014), which has been compiled in the framework of the Global Earthquake Model (GEM, 2015) and includes all historical, early and modern instrumental events up to 2006 with unified moment magnitude.For the Makran region, the authors suggest a magnitude of completeness of M C = 5.5 after 1920 and M C = 4.5 after 1963.We select only those events (colored circles on left panel of Fig. 1, listed in Supplement 7) lying in the area covered by the plate interface, which clearly excludes the non-subduction seismicity to the west and east (Fig. 1).
The seismicity is determined using the procedure by Kijko and Smit (2012) using their Matlab code "aue.m" and applying magnitudes of completeness following Zare et al. (2014) using the events from 1919 to 1963 and from 1964 to 2006 with an assumed magnitude error of 0.2.The resulting b value is 0.82.The code "aue.m" also generates cumulative seismicity rates (lambda) based on Gutenberg-Richter-Bayes (GRB) statistics which include uncertainty for the magnitude-frequency relations (right panel of Fig. 1).We use this distribution for the generation of the reference synthetic catalog (RSC).An observational time span of less than 100 years is actually not sufficient to estimate the largest possible earthquake, with a return period of possibly more than 1000 years.This is why we systematically assess the effect of maximum magnitude assumption on tsunami hazard.
The plate interface geometry we based our seismicity analysis and the synthetic catalog on is described and shown in Appendix A. It clearly excludes the non-subduction-related seismicity to the west and east (Fig. 1).

Synthetic catalog
The synthetic catalogs each span 300 000 years and are created by randomly generating events within the area spanned by the plate interface geometry with magnitudes according to the inferred magnitude-frequency relation (and exponential inter-event time distribution).All the events are assumed to be thrust events at the subduction interface.Of course, this assumption would be unacceptable for a seismic hazard analysis, since part of the seismicity is related to crustal and intraplate earthquakes in the slab.However, in the case of a tsunami hazard assessment it is not such an issue: all the events on land do not cause any tsunami at all, whereas most of the large, tsunamigenic events in the sea can reasonably be assumed to be thrust interplate events (strike-slip events have lower impact).The reason why we do include the area on land is to get a better picture of the seismicity (which is low and furthermore not well recorded) for the determination of the magnitude-frequency relation.The maximum magnitude for the RSC is set to 9.0.It is important to stress that this value is not the most likely one, but it has to be considered as end-member assumption.
In addition to the RSC, other catalogs are generated by varying the Gutenberg-Richter a and b values, a combination thereof resulting in a rotation about the rate at M = 6, and, most importantly, the maximum magnitude (Fig. 2).
The hypocenters of the catalog events are projected to the subduction interface, which consists of subfaults of about 20 km edge length.These are activated based on empirical scaling relations for earthquake extension and slip distribution.More details of the algorithm and examples are given in Appendix C.

Tsunami model
For every subfault, the sea floor deformation resulting from 1 m dip slip (rake angle is set to 90 • ) is calculated for the homogeneous elastic half-space and used as input for the tsunami code easyWave (Babeyko, 2012), which propagates the initial condition to compute wave height time series for points along the coast in about 50 m water depth (tsunami Green's functions).Up to this water depth it is possible to assume linear conditions for our purpose (Miranda et al., 2014),  , 2014).The brown lines are the profile and deformation front by Smith et al. (2013), in grey the subduction interface.The color of the selected events corresponds to year, and size is by scaling relations (Blaser et al., 2010).Right: inferred magnitude-frequency relation.Solid: Gutenberg-Richter-Bayes (used), grey: Gutenberg-Richter (for comparison).so that the unit time series can be combined according to the slip distribution.This enables us to compute a large number of scenarios with little computational effort.The maximum of each scenario time series is extrapolated to the shore (1 m water depth) using Green's law (Kamigaichi, 2009) to obtain peak coastal tsunami amplitude (PCTA) at points of interest (POI).Bathymetry is from GEBCO (2008) with 30 arcsec resolution.

Tsunami hazard
For each point of interest along the coast (or collection of points) the number of exceedances of a certain wave height level is counted.The percentiles are calculated using the inverse cumulative Poisson distribution.For instance, if there are 150 exceedances at some point for some wave height, the 15th and 85th percentiles are 137 and 163 (Matlab: poissinv([0.15,0.85], 150)), which means that 70 % of random realizations of a synthetic catalog with identical parameterization will fall within these numbers.We test this in Appendix B.
Dividing the number of exceedances by catalog time span yields the annual rate R. The probability of exceedance (POE) for N years is computed using Poisson statistics: The probabilistic tsunami height (PTH) for time T is defined as the tsunami height which is on average exceeded once in T years (return period) (Bernard and Robinson, 2009;Sørensen et al., 2012).It is computed by interpolating wave height as a function of 1/R at T .Following the above formula, the PTH is exceeded with a probability of 1−exp(−1) ≈ 0.63 in T years one or more times.
The hazard for the cities is not calculated based on a single calculation point, but the nine closest points, resulting in somewhat higher numbers.one is much higher, about 21 and 9 m, respectively (circles with bars).The central panel shows PTH as a function of time for the whole coastline and for some cities along the coast.Jiwani (at 61.7 • Lon) has the highest hazard, since it is exposed to more events than the cities of Jask or Ormara at the boundaries of the subduction zone.The circles correspond to the circles in the left panel.The right panel shows the POE at one or more locations as function of PCTA.The squares denote POE for 2 m/50 yr, 8 m/500 yr and 15 m/5000 yr and are used for comparison in the summarizing Figs.5-6.The dashed lines are the 15th and 85th percentiles and reflect fluctuations due to stochastical catalog generation.
Figure 4 shows the same plots for the exposed coastline of Oman.The hazard for Muscat (23.6 • Lat) and Sur (22.6 • Lat) is comparable to the one for Chabahar in Iran, but the cumulative hazard for Oman is lower, since only the events in the western part have a heavy impact on Oman.
The effect of assumed maximum magnitude for the synthetic seismic catalog on tsunami hazard is shown in Fig. 5, for Iran and Pakistan on the left side and for Oman on the right side, as PTH in the upper and as POE in the lower panels.The low-level short-term hazard is dominated by the large number of smaller events, while the high-level long-term hazard very strongly depends on the maximum magnitude.For instance, for Iran and Pakistan, the probability of exceeding 15 m in 5000 years increases from 21 % for M Max = 8.2 to 75 % for M Max = 8.6 and 93 % for M Max = 9.0.
Similarly, Fig. 6 shows the dependence of tsunami hazard on magnitude-frequency relations (for Iran and Pakistan and Oman together).Throughout Figs. 3 to 6, circles and squares correspond.Numerical values are listed in Supplements 1-4.

Discussion
Due to the rather low seismic activity, the tsunami hazard for Makran is generally lower than for other, seismically more active subduction zones, especially concerning short-term hazard.However, on the long-term, events with potentially catastrophic impact lead to significant hazard, meaning that risk mitigation efforts and early warning measures should be further developed and implemented.Our estimations of tsunami hazard should not be considered to be final.We discuss some of the factors affecting our results.
Stability of the results from the stochastically generated synthetic seismic catalogs -It is estimated by computing the 15th and 85th percentiles and is high due to the length of the catalogs of 300 000 years.A more detailed discussion including additional test catalogs is given in Appendix B and Supplement 5. Earthquake interaction in space and time -The temporal and spatial interplay of earthquakes is a complicated process which has to be considered, for instance, if simulations of the earthquake cycle are carried out.In this study however, we can assume that for a long time span, the average slip will be similar across the interface (homogeneous plate coupling assumed).Temporal clustering is actually irrelevant for our hazard computations, since only the total time span of the catalog is used therefore.
Assumptions for scaling relations and slip distribution -These could have high impact on single events, where more complicated things like asperities can arise.Averaged over the whole catalog, the effect should be moderate.The effect of modified scaling relations and slip distribution is shown in Appendix C and Supplement 6.
Projection of all events to the plate interface -Although this procedure would be inadequate for seismic hazard assessment, for tsunami hazard assessment it is legitimate and leads to a more conservative view as discussed in Sect.2.2.
Splay faulting -We did not explicitly include geometric splay faults, which could increase the hazard.However, their effect is partly reproduced by our slip-distribution algorithm, as explained in Appendix C.
Tsunami modeling -We did not compute inundation maps and run-up, something which would be useful for the analysis of a few specific hypothetic tsunami events.Instead, we estimated coastal wave heights as described in Sect.2.3.Using Green's law for the estimation of tsunami coastal impact has become common practice when computing wave propagation for a very large number of scenarios.The reasoning is as follows.The number of grid elements is not able to resolve the shoaling of the wave when it approaches the coast and the wave height at the coast will be underestimated.The use of nested grids or non-regular finite-element meshes dramatically increases the computational cost.Instead, Kamigaichi (2009) demonstrated that Green's law projection from offshore positions gives coastal wave heights similar to those computed on a much finer grid.This approach is also implemented, for example, in the Japanese early warning system.
Submarine landslides -These could possibly be triggered by earthquakes and increase the hazard.Technically, they could be handled in modeling, but due to even lower rate of occurrence than for earthquakes it is difficult to assign probabilities.
In our opinion, the largest influence stems from uncertainty of the maximum magnitude, which we tried to assess in this study.Actually, to our knowledge, this is the first PTHA sensitivity analysis concerning maximum magnitude.
We now compare our results to findings by Heidarzadeh and Kijko (2010).For Jiwani, they report probabilities of exceeding 1, 3 and 5 m over a 50 year period of about 44, 18 and 18 %.The authors refer to their publication as a "first generation" study, since it contains only three earthquakes with identical magnitudes of M W = 8.1 (our catalogs contain about 20 000 events).This rough discretization leads to the same value for 3 and 5 m exceedance probabilities.Another limitation of the above study is that the PTHA is only valid for times that are not too far away from the typical recurrence time of a M W = 8.1 event, while we cover the full spectrum of magnitudes and thus hazard times.That all three events affect Jiwani also implies that these values should not be compared to our results for a specific location, but to our probability of exceedance anywhere within some coastal region.Our values for specific locations are clearly smaller (see central panels of Figs. 3 and 4).For a maximum magnitude of M Max = 8.2 our results are 46, 24 and 14 %, for the coastline of Iran and Pakistan, which is in good agreement with the result above.For M Max = 8.6 the values increase rather moderately to 49, 28 and 19 %.As mentioned in Sect.3, this short-term hazard is dominated by smaller events, so that higher M Max has limited effect, but the long-term hazard is heavily affected by M Max .
An interesting question concerns the possible segmentation of the Makran subduction zone.The seismicity varies strongly across the Makran region.While large thrust events are known to have occurred in the eastern part, there is no proof of historic events in western Makran (Musson, 2009).Byrne et al. (1992) discussed whether this should be interpreted as aseismic sliding in the west, or rather complete locking of the interface (with opposite hazard consequences).They concluded that the occurrence of uplifted Holocene marine terraces along the coast as well as the continuity of the deformation front indicate that the latter is more probable.However, even assuming locking, it is still not clear whether a giant event rupturing the whole plate interface at once is possible or whether the segmentation would favor sub-events involving only part of the interface.In this study we did not consider segmentation of the subduction zone; further studies with refined seismicity models taking segmentation into account would be desireable.
Looking at the POE for 5000 years in the lower left panel of Fig. 5, we see that there is a very strong dependence on M Max , which could pose a huge problem for decision makers.However, we have to stress that not all M Max are equiprobable.Since we know that an event with M W = 8.1 has already occurred, the values below M Max = 8.1 can be neglected, and are shown only to illustrate the general principle, which is also the reason why we show the values to the far right side.An event with M W = 9.4 cannot be reasonably fit into a subduction zone of that size, and also M W = 9.2 is only feasible under very extreme and unrealistic assumptions.Nevertheless, even in the realistic range, the POE for 15 m in 5000 years still varies between maybe 40-80 %.This reflects the steep transition necessarily occurring for the probability of something which is almost impossible (e.g., exceeding 15 m with magnitudes below M Max = 8) to something that is probable (and also almost certainly happens in the long run for large M Max ).For this reason it is maybe more useful to look at PTH instead of POE, which is more linearly dependent on M Max .For the corresponding range, for 5000 years it varies from around 14-18 m.This variation of 25 % should be manageable for decision makers.For smaller hazard times, both PTH and POE are less dependent on M Max .
We close by listing some general recommendations beyond hazard analysis on how to increase resilience against tsunami as education of the population, planning of evacuation procedures, protection of high-risk facilities and installation of efficient early-warning systems.The traditional approach for tsunami early warning is based on seismic methods, but especially for Iran and Pakistan, a GNSS (Global Navigation Satellite System) based approach should be considered (Sobolev et al., 2007;Hoechner et al., 2008Hoechner et al., , 2013;;Babeyko et al., 2010).

Appendix A: Subduction geometry
The Makran geometry is neither available via the RUM model (Gudmundsson and Sambridge, 1998) nor Slab1.0 (Hayes et al., 2012).We used a profile of the slab by Smith et al. (2013, Fig. 2a), which is based on a combination of seismic reflection lines, and which the authors employed for their thermomechanical modeling.The profile is shown in map and side view in Fig. A1 (solid red line).We slightly shifted the profile (dashed line) and rotated it around 61.75 • E/38.0 • N so as to match the deformation front (blue line) from the same publication (Fig. 1 therein).This results in the 3-D geometry as shown in Fig. A1 in green (rectangles are the earthquake-and tsunami subfaults for the Green's functions).The lower edge is at 70 km depth.This approach by no means provides an accurate, final model of the subduction interface.Many more data and studies are highly desirable to get a better understanding of the structure (as well as the activity) at the Makran.Hazard models based on other geometry assumptions could be made, but taking into account geometric complications or variations is outside the scope of the present study, where the focus is on the consequences of maximum magnitude assumptions.In addition to Fig. 2, where we plotted the epicenters of the first 5000 events of the synthetic catalog, in Fig. B1, we plot the epicenters for the full 300 kyr for magnitude above M W = 8 and M W = 8.5, showing that indeed the coverage of the area is fairly complete, especially considering the extension of the large events.In order to further assess the stability of the results with respect to the stochastic nature of the synthetic catalog, we generated four additional test catalogs with the same parameterization as the RSC.The results based on these catalogs are listed in Supplement 5 and shown in Fig. B2.73 % of the results fall within the 15th and 85th percentiles (70 % expected).

Appendix C: Scaling relations and slip distribution
The most frequently used empirical relation for earthquake extension is probably that by Wells and Coppersmith (1994).However, we use the one by Blaser et al. (2010) (reverse orthogonal), because it is based on a more recent data set and specifically considers subduction events.
In our algorithm, a subfault is activated if the dip and strike components of the distance of its center to the hypocenter fall within the empirical scaling relations.Along the dip, the strength of the activation is determined by the formula by Geist and Dmowska (1999) (symmetrical), which is analytically derived based on physical assumptions and leads to a bell-shaped slip distribution.Within 20 % of the lateral edges (along strike) we linearly taper the slip to zero.This means that in general, the randomly selected hypocenter is (almost) the barycenter of the slip distribution.If, however, the hypocenter falls close to an edge of the plate geometry, the missing moment due to fewer activated subfaults is compensated by higher slip on the other subfaults, and the slip is not tapered to zero at the outer edge.This can partly mimic the effect of splay faults at the surface.In this case, the event does not strictly follow the scaling relations and the hypocenter does not correspond exactly to the barycenter.Some examples for events of varying magnitude are shown in Fig. C1.In order to assess the effect of the empirical scaling relations and slip distribution on tsunami hazard, we generated additional catalogs based on the reference catalog.For each event we modified length and width randomly up to 50 %, but so as not to exceed 50 % variation in their product and quotient (Fig. C2).In an other version we perturbed the slip distribution so as to produce irregular slip.This was done by multiplying the slip distribution with a spatially low-pass filtered uniform random distribution (Fig. C3).Finally we applied both modifications together.A comparison of the reference results and the modified versions is given in Supplement 6.
A. Hoechner et

Figure 2 .
Figure2.Left: first 5000 events in the reference synthetic catalog color-coded by magnitude (larger on top).The whole catalog spanning 300 000 years contains about 20 000 events.Right: parameter range tested for the magnitude-frequency relation.

Figure 3 .Figure 4 .
Figure 3. Left: probabilistic tsunami height (PTH) for Iran and Pakistan along Longitude (colored lines) and for the whole coast at once (circles with 15th and 85th percentile bars).Center: PTH as function of time for selected cities (color) and the coastline as whole (black).Right: probability of exceedance (POE) as function of peak coastal tsunami amplitude (PCTA) for different time periods.Dashed: 15th/85th percentiles.

Figure 6 .
Figure 6.Influence on PTH (upper panels) and POE (lower panels) from varying a parameter (seismic rate) (left) and b parameter (right) around M = 6.

Figure A1 .
Figure A1.Profile (red) and deformation front (blue) as from Smith et al. (2013) and constructed 3d geometry (green) as used in the present study.

A
. Hoechner et al.: Probabilistic tsunami hazard assessment Appendix B: Synthetic catalog stability

Figure B1 .Figure B2 .
Figure B1.Same plot as in Fig.2, but for the whole 300 kyr of the reference synthetic catalog for magnitudes above 8 and 8.5.

Figure C1 .Figure C2 .
Figure C1.Illustration of the empirical scaling relations and slip distribution.
al.: Probabilistic tsunami hazard assessment Examples for the randomization of the slip distribution (and rake angle).The Supplement related to this article is available online at doi:10.5194/nhess-16-1339-2016-supplement.