Probabilistic tsunami hazard assessment for the Makran subduction zone using stochastic slip models

Abstract

The Makran subduction zone has generated damaging tsunamis, most notably the 27 November 1945 event, but its future hazard remains difficult to quantify because plate coupling, segmentation, splay faulting, shallow slip amplification, and recurrence rates are poorly constrained. We present a probabilistic tsunami hazard assessment for the Makran coasts of Iran, Pakistan, Oman, and western India using stochastic megathrust slip models embedded in a logic-tree framework. Earthquake sources from Mw 7.8 to 9.2 were sampled across western, eastern, and full-margin rupture geometries. Slip fields were generated as spatially correlated random fields conditioned on seismic moment, asperity depth, shallow-slip amplification, and long-term moment balance. Seafloor deformation was computed with elastic dislocation solutions and propagated with a depth-averaged nonlinear shallow-water solver over nested bathymetric grids. The ensemble contains 48,000 earthquake scenarios and 2,400 landslide-screening scenarios; hazard curves were produced for 32 coastal sites and for gridded offshore points at 20 m depth. The median 475-year maximum tsunami amplitude is 1.1-2.4 m along exposed Makran coasts and 0.25-0.75 m along the Gulf of Oman, while the 2,475-year amplitude reaches 3.8-7.1 m at selected headlands. Uncertainty is dominated by recurrence-rate branching at long return periods and by shallow slip amplification at near-source sites. Validation against published 1945 tsunami interpretations favours eastern-segment rupture with heterogeneous shallow slip but does not rule out multi-segment events. The study supports risk-informed planning based on hazard curves rather than a single worst-case scenario, while emphasising that inundation maps require local topography, tide, roughness, and exposure data beyond the offshore hazard estimates reported here.

Introduction

The Makran subduction zone is one of the least instrumented but most consequential tsunami sources around the northwestern Indian Ocean. The Arabian plate subducts beneath Iran and Pakistan, and the margin has produced large thrust earthquakes, aseismic slip, coastal deformation, and the damaging 1945 tsunami [7,8,9,10]. The historical record is sparse compared with Japan, Chile, or Sumatra, so hazard estimates depend strongly on how geological and geodetic uncertainty is translated into source models.

Previous Makran tsunami studies have moved from deterministic scenarios to probabilistic assessments and stochastic rupture models [1,2,3,4,5,6,19]. That shift is important. A single Mw 8.3 or Mw 9.0 scenario can support emergency planning, but it cannot answer questions about annual exceedance probability, return-period amplitude, or which uncertainty branch controls design levels. Probabilistic tsunami hazard analysis provides the needed framework by combining source recurrence, rupture variability, propagation, and exceedance statistics [11,12].

The main unresolved issue is shallow slip. Tsunami height is controlled not only by earthquake magnitude, but also by where slip occurs relative to the trench, whether splay faults participate, and how heterogeneous asperities focus uplift. Stochastic slip models capture some of this variability through correlated random fields and moment constraints [13,14,15]. However, they can also create physically unrealistic scenarios if slip, rupture area, and recurrence are sampled independently. Recent work therefore emphasises long-term moment balance and explicit uncertainty in shallow slip amplification [18].

This article reports a fictional but realistic PTHA workflow for the Makran margin. It combines source branches from Makran-specific literature, stochastic slip ensembles, elastic seafloor deformation, nonlinear shallow-water propagation, and site hazard integration. The output is an offshore hazard product at 20 m water depth, not a final inundation map. That distinction matters because nearshore amplification and overland flow depend on local bathymetry, coastal topography, tide, and roughness.

Tectonic and historical setting

The Makran margin is commonly divided into western and eastern segments. The eastern segment ruptured in 1945, producing a tsunami that affected the coasts of present-day Pakistan, Iran, Oman, and India. Source studies of the 1945 event infer Mw about 8.1-8.3, with heterogeneous slip and coastal deformation compatible with a shallow megathrust source [6,7]. Numerical modelling of the event has reproduced broad arrival-time and amplitude patterns but also shows strong sensitivity to the assumed rupture geometry [8].

The western segment is more ambiguous. Historical seismicity is lower, but low seismicity does not necessarily imply low hazard. Byrne, Sykes, and Davis argued that aseismic slip and great thrust earthquakes must both be considered along the plate boundary [10]. More recent seismotectonic reviews emphasise uncertain segmentation, sparse offshore geodesy, sediment thickness, and possible splay faulting [9]. These uncertainties are exactly the kind of problem a logic tree is meant to expose rather than hide.

Makran hazard studies have reached different absolute amplitudes because they use different recurrence models, magnitude limits, and rupture geometries. The 2011 northwestern Indian Ocean PTHA, later comprehensive Makran assessments, and recent stochastic-source studies all agree on one qualitative point: near-source coasts can experience damaging waves from plausible megathrust ruptures, and the tails of the hazard curve depend on rare multi-segment or high-shallow-slip events [1,2,3,4,5].

Our source model therefore avoids treating the 1945 event as the only template. It uses the 1945 data as a validation anchor and samples broader rupture families for future hazard.

Source-zone logic tree

The logic tree contains four main source-zone branches: eastern Makran, western Makran, full-margin rupture, and splay-assisted rupture. The eastern and western branches represent ruptures confined mainly to one segment. The full-margin branch represents rare events crossing the central structural transition. The splay-assisted branch adds shallow deformation on an imbricate fault near the accretionary wedge, motivated by previous Makran studies that evaluate splay fault contributions [2].

Magnitude-frequency rates were assigned using three recurrence branches. The low-rate branch follows a conservative historical-seismicity interpretation; the central branch balances long-term plate convergence with partial seismic coupling; and the high-rate branch permits larger locked fractions in the eastern and western segments. Branch weights were chosen to be transparent rather than optimised to one preferred hazard result: 0.25, 0.50, and 0.25 for low, central, and high recurrence. The full-margin branch carries small annual probability but high consequence. The recurrence treatment was also checked against broader Indian Ocean tsunami recurrence estimates so that Makran rates were not interpreted in isolation [20].

Maximum magnitude was varied between Mw 8.6 and 9.2 depending on source-zone branch. We did not allow Mw 9.2 events on the isolated eastern branch because the rupture area would exceed the segment width used in the structural model. Conversely, the full-margin branch permits Mw 9.1-9.2 events only under high coupling and large down-dip width. This coupling between magnitude and geometry prevents the ensemble from generating impossible ruptures.

The logic-tree design follows the broader PTHA principle that epistemic uncertainty should be represented by alternative branches, while aleatory rupture variability should be represented by scenario sampling within a branch [11,12].

Stochastic slip generation

For each rupture, we generated slip as a two-dimensional lognormal random field on the fault plane. Correlation lengths scale with rupture dimensions, and the mean slip is set by seismic moment. The random-field approach follows the general earthquake-slip framework of Mai and Beroza [14], modified for tsunami applications by enforcing non-negative slip, smoothing unrealistic cell-to-cell jumps, and conditioning the shallow slip branch.

Three shallow-slip branches were used. The neutral branch applies no systematic depth trend. The amplification branch increases mean slip above 15 km depth by a factor sampled between 1.4 and 2.6 while conserving moment by reducing deeper slip. The suppressed-shallow branch moves asperities down-dip, representing ruptures whose shallow interface is weakly coupled or decoupled. These choices reflect the strong influence of shallow slip on tsunami hazard in subduction zones [18].

Each magnitude-source-zone branch used 400-900 stochastic slip realisations, giving 48,000 earthquake scenarios in total after recurrence weighting. Slip fields were rejected if their peak slip exceeded five times branch mean slip, if the rupture centroid fell outside the intended segment, or if moment conservation failed by more than 1%. The rejection rate was 6.8%, mostly from high-shallow-slip full-margin scenarios.

The 1945 calibration subset used the same slip generator but was conditioned to Mw 8.1-8.3 and eastern-segment rupture. This allowed us to compare the model family with published 1945 waveform and deformation interpretations without tuning the entire hazard model to that event [6,7].

Seafloor deformation and propagation

Static seafloor deformation was computed from rectangular subfault dislocations in an elastic half-space using the Okada formulation [16]. The older Mansinha-Smylie displacement solution was used as a cross-check for selected scenarios [17]. Horizontal displacement contributions on sloping bathymetry were included for the main runs because the Makran accretionary wedge has steep local gradients and thick sediments.

Initial sea-surface displacement was assumed equal to coseismic seafloor displacement after applying a 90 s rise-time filter. We did not model dynamic rupture or dispersive wave generation. This is acceptable for the offshore hazard metrics used here, but it may underestimate short-wavelength components near steep bathymetry or landslide sources.

Tsunami propagation used a finite-volume nonlinear shallow-water solver on nested grids with 30 arcsec outer resolution and 3 arcsec nearshore resolution around selected coastal sites. The numerical configuration was benchmarked against standard tsunami-current cases using the same class of shallow-water methods evaluated in Tsunami-HySEA benchmarking studies [21,22]. Wetting and drying were disabled for the primary hazard curves because the output is extracted at 20 m depth offshore rather than on land.

Gauge points were placed offshore of Chabahar, Gwadar, Pasni, Ormara, Karachi, Muscat, Sur, Khasab, and selected smaller headlands. For each scenario we stored maximum positive wave amplitude, first-arrival time, and the time-integrated squared velocity proxy. Maximum amplitude is the primary hazard measure in this article; velocity proxies are provided in the supplementary data for future harbour and current studies.

Hazard integration

Annual exceedance rates were computed by summing scenario probabilities for amplitudes above each threshold. For each coastal point, the hazard curve gives lambda(A > a), where A is maximum positive offshore amplitude and a is the amplitude threshold. Return-period amplitudes are then obtained by inverting the annual exceedance curve. We report 475-year, 975-year, and 2,475-year amplitudes because these are familiar engineering intervals, not because tsunami recurrence is periodic.

Epistemic uncertainty is reported as fractile curves across logic-tree branches. Aleatory rupture variability is represented within each branch by stochastic slip realisations. This separation follows standard PTHA practice [11,12]. We also compute a deaggregation table showing which magnitude, segment, and shallow-slip branch contributes most to each site and return period.

Scenario weights were normalised so that each recurrence branch satisfies its moment-rate constraint over a 20,000-year synthetic catalogue. The catalogue length is only a numerical integration device; it should not be read as a literal simulation of future history. The final hazard curves use annual rates, not catalogue counts.

A landslide-screening calculation was included but not merged into the main hazard curve. The Makran accretionary wedge contains thick sediment and potential slope failures, yet landslide recurrence and volume distributions are much more uncertain than megathrust recurrence. We therefore report landslide-sensitive sites qualitatively and leave full multi-source integration for future work.

Validation against the 1945 event

The 1945 event is the only strong regional calibration point, but the observations are incomplete. We compared the conditioned eastern-segment subset with published waveform and coastal deformation interpretations [6,7] and with numerical modelling of the 1945 tsunami [8]. The validation used arrival-time windows and broad amplitude classes rather than exact point amplitudes, because historical run-up observations include large local and documentary uncertainty.

The best-fitting stochastic realisations had Mw 8.2, rupture length about 160-210 km, peak slip of 6-9 m, and shallow slip concentrated offshore of Pasni and Ormara. Uniform-slip models reproduced first arrival times but under-predicted site-to-site amplitude variability. This supports the use of heterogeneous slip for hazard analysis, consistent with stochastic 1945 hazard studies [6].

The validation does not prove that future eastern Makran ruptures will look like 1945. It only checks that the source generator can produce known-event-like wavefields without special tuning. Several plausible 1945 realisations also produced much smaller waves in Oman than some historical accounts suggest, indicating that local bathymetry, tide, and reporting uncertainty remain important.

Because the western Makran has no equivalent modern tsunami calibration event, its hazard curves remain more dependent on tectonic and recurrence branches. We keep this asymmetry visible in the uncertainty bands rather than forcing western and eastern segments to share the same confidence.

Offshore hazard results

The median 475-year offshore amplitude at 20 m depth is highest along the exposed Pakistan and southeastern Iran coasts. Pasni, Gwadar, and headlands west of Chabahar have median 475-year amplitudes between 1.6 and 2.4 m. Karachi has lower median amplitude, 0.55 m, because of distance and shelf geometry, but the 95th percentile remains above 1.5 m for high-shallow-slip eastern ruptures.

The Gulf of Oman sites have lower but non-negligible hazard. Muscat and Sur have median 475-year amplitudes of 0.25-0.75 m and 2,475-year amplitudes of 1.2-2.6 m. The largest Oman amplitudes come from full-margin and western-segment ruptures rather than from 1945-like eastern events. This agrees qualitatively with previous Makran hazard assessments that identify strong spatial variation across the northwestern Indian Ocean [1,3,5].

At 2,475-year return period, the median hazard becomes strongly controlled by full-margin and shallow-amplified branches. Exposed Makran headlands reach 3.8-7.1 m offshore amplitude in the central branch and exceed 8 m in the 95th percentile for several sites. These numbers should not be interpreted as onshore run-up. Local inundation could be larger or smaller depending on bathymetric focusing, harbour resonance, tide, and coastal topography.

Arrival times are short for near-source sites. Median first-arrival times are 18-27 min for Pasni and Gwadar, 35-52 min for Chabahar, 70-105 min for Muscat and Sur, and longer than 90 min for Karachi in most scenarios. The short near-source warning time reinforces the importance of pre-event planning rather than response based only on far-field confirmation.

Hazard deaggregation

Hazard deaggregation shows that different sites are controlled by different parts of the source model. For Gwadar at the 475-year level, eastern-segment Mw 8.2-8.6 ruptures contribute 61% of exceedance rate, with shallow-amplified slip contributing just under half of that. At the 2,475-year level, full-margin Mw 8.8-9.1 ruptures become the largest contributor.

For Chabahar, western-segment and full-margin branches contribute more strongly than for Gwadar. The site is therefore more sensitive to the poorly constrained western segment. For Muscat, western and full-margin sources dominate at all return periods above about 500 years. This explains why a hazard model calibrated only to 1945 can understate Omani tail hazard.

The deaggregation also shows the influence of stochastic slip within the same magnitude. At Pasni, the 90th percentile of amplitude for Mw 8.4 eastern ruptures is 2.3 times the 10th percentile. Magnitude is therefore a poor single proxy for tsunami severity. Slip centroid, shallow amplification, and along-strike asperity location all matter.

These results are consistent with the PTHA literature: hazard is not a deterministic function of maximum magnitude, and multiple source branches can dominate different sites or return periods [11,12,18].

Sensitivity analysis

The largest epistemic sensitivity is recurrence rate. Moving from the central to high recurrence branch increases 2,475-year amplitudes by 18-42% depending on site. The effect is strongest for western Makran and Oman sites, where the historical record provides the weakest constraint. This is unsurprising but important: better offshore geodesy and palaeotsunami constraints would directly reduce hazard uncertainty.

The second sensitivity is shallow slip amplification. At near-source Pakistan sites, activating the shallow-amplification branch raises the median 475-year amplitude by 22-38% and the 2,475-year amplitude by up to 51%. At farther-field sites, shallow slip still matters but is partly averaged by propagation distance. Recent subduction-zone PTHA work shows why this uncertainty should be represented explicitly rather than absorbed into a generic rupture model [18].

Fault dip and down-dip width have site-dependent effects. Steeper dip increases near-trench uplift but can reduce broad uplift area. Wider down-dip ruptures increase moment for the same slip but may shift uplift landward, changing offshore amplitudes. The sensitivity is smaller than recurrence and shallow slip for most sites, but it affects the western-segment tail.

Numerical-grid sensitivity was below 8% for offshore gauges outside complex harbours. Near headlands, differences between 6 arcsec and 3 arcsec nearshore grids reached 15%. This is why we stop at offshore hazard. Inundation mapping should use higher-resolution topobathymetry and local calibration.

Discussion

The main practical result is that the Makran tsunami hazard cannot be summarised by one 1945-like scenario. The 1945 event is essential for validation, but the hazard at long return periods is shaped by rare full-margin ruptures, western-segment uncertainty, and shallow slip. This agrees with the direction of recent Makran stochastic PTHA work [1,2,19].

For emergency planning, the offshore hazard curves support two complementary products. The first is a frequent-to-moderate event set for warning exercises and evacuation timing. The second is a low-probability high-consequence set for critical infrastructure and land-use screening. Mixing these into one "maximum credible tsunami" can obscure probability and encourage brittle design assumptions.

The results also show where research would most improve the hazard model. Offshore geodetic constraints on coupling, marine seismic imaging of splay faults, palaeotsunami records along Oman and Iran, and better sediment-slope stability data would all reduce uncertainty. More numerical propagation alone will not solve the source problem.

Finally, the modelling choices are deliberately transparent. We do not claim that stochastic slip fields are the truth. They are a disciplined way to sample rupture complexity when observations are sparse. The value of the framework is that assumptions can be changed and their effect on hazard curves can be measured.

Limitations

The first limitation is source uncertainty. The western Makran recurrence rate, coupling fraction, and segmentation are not well constrained. We represent this uncertainty in the logic tree, but a broad logic tree is not a substitute for data. The western-tail hazard should be read as plausible rather than tightly estimated.

The second limitation is the offshore metric. Maximum amplitude at 20 m depth does not equal run-up, flow depth, or damage. Coastal amplification, tide, harbour resonance, wave breaking, roughness, and building exposure can change the risk substantially. Site-specific inundation studies must follow before design elevations are set.

The third limitation is source physics. The simulations use static seafloor deformation and kinematic rise-time smoothing. Dynamic rupture, poroelastic sediment response, landslide coupling, and dispersive waves are not fully represented. These effects may matter for shallow wedge deformation and near-field short-period waves.

Finally, landslide tsunamis are screened but not integrated probabilistically. The Makran wedge is sediment rich, and landslide sources may be important locally. Their recurrence and volume distributions remain too uncertain for the same treatment as megathrust earthquakes in this article.

Conclusion

A stochastic-slip probabilistic tsunami hazard assessment for the Makran subduction zone indicates substantial near-source offshore hazard along the Pakistan and southeastern Iran coasts and lower but significant hazard along the Gulf of Oman. The 1945 tsunami is best reproduced by heterogeneous eastern-segment shallow slip, but long-return-period hazard is controlled by full-margin, western-segment, and shallow-amplified scenarios. Recurrence uncertainty and shallow-slip amplification dominate the epistemic spread.

The article supports hazard communication through site-specific exceedance curves, deaggregation, and scenario families rather than through a single worst-case map. For practical risk reduction, the next step is to couple these offshore hazard curves to high-resolution inundation modelling and exposure analysis for priority ports, coastal towns, and critical infrastructure corridors.

Data and code availability

The supplementary archive contains source-zone polygons, branch weights, stochastic slip realisations, tsunami input files, offshore gauge coordinates, hazard curves, deaggregation tables, and plotting scripts. Bathymetry files are referenced by source and checksum rather than redistributed where licensing restricts direct packaging.